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ABSTRACT 


In  this  study,  the  time  temperature  profile  of  a  missile  exposed  to  fire  in  a 
compartment  adjacent  to  the  missile  magazine  is  examined.  The  study  required  the 
development  of  a  heat  transfer  model  based  on  the  geometry  and  thermophysical 
properties  of  a  new  concept  for  a  vertical  launching  system,  the  Concentric  Canister 
Launcher  (CCL).  Different  fire  scenarios  are  analyzed  by  the  model  to  predict  the  time  it 
takes  to  reach  a  critical  value  or  “cook-off’  temperature  of  the  missile’s  propellant  and 
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I.  INTRODUCTION 


A.  BACKGROUND 

The  history  of  naval  combat  has  shown  that  there  is  a  demonstrated  need  to 
understand  how  shipboard  fires,  and  their  containment,  affect  ship  survivability.  Credible 
fire  threat  models  are  particularly  helpful  in  assessing  the  adequacy  of  fire  protection 
systems,  and  in  providing  baseline  data  to  make  sound  damage  control  decisions  to 
combat  conflagrations  which  may  jeopardize  munition  magazines.  Fires  initiated  by  the 
penetration  and  detonation  of  Excocet  missiles,  and  sustained  with  their  unspent  missile 
propellent,  were  experienced  by  the  USS  STARK  (FFG  3 1)  in  the  Arabian  Gulf  and 
British  ships  during  the  Falkland  conflict.  In  all  instances,  ship’s  survivability  was 
seriously  jeopardized.  In  the  Falklands,  a  British  ship  was  abandoned  and  later  lost.  In 
each  instance,  munitions  magazines  were  threatened  by  fire  that  spread  into  adjacent 
compartments.  Drastic  measures  were  required  to  cool  the  munitions,  and  prevent  the 
detonation  of  high  explosive  material  and  the  ignition  of  solid  propellent.  [Ref.  1] 

B.  CCL  CONCEPT 

Naval  Surface  Warfare  Center  Dahlgren  Division  (NSWC)  [Ref.  2]  is  currently 
examining  a  new  concept  for  a  future  vertical  launching  system,  the  Concentric  Canister 
Launcher  (CCL).  This  new  launcher  will  provide  greater  firepower,  flexibility,  and  lower 
cost  based  on  its  modular  design,  than  the  current  MK  41  VLS.  The  CCL  was  originally 
proposed  to  replace  the  MK  26  Guided  Missile  Launching  System  (GMLS)  aboard  CGN 
36,  DDG  993,  CG  47-51,  and  non  VLS,  DD  963  classes  to  enhance  capability  with 
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increased  magazine  capacity.  Additionally,  the  CCL  would  be  adaptable  to  forward  fit 
new  ships  such  as  future  DDG-51's,  SC-21,  Arsenal  and  Swath  hull  designs. 

A  schematic  of  the  CCL  is  given  in  Figure  1 .  It  consists  of  two  concentric 
cylinders  joined  by  sets  of  dual  longerons.  One  end  is  open  and  the  other  is  sealed  with  a 
hemispherical  cap.  The  annulus  between  the  inner  and  outer  cylinder  provides  gas 
management  by  ducting  discharge  exhaust  combustion  gases  at  launch.  Each  canister 
carries  one  missile  and  fits  into  a  ship’s  weapon  module.  The  weapons  module  consists  of 
a  test  tube  type  rack  in  which  each  missile  canister  is  held  by  a  shock  collar  at  the  main 
deck.  Schematics  of  the  ship’s  weapon  module  and  shock  collar  are  given  in  Figures  2 
and  3,  respectively.  The  potential  advantages  of  this  concept  over  the  current  MK  41  are: 

•  Round  pressure  vessels  and  one  shot  design  allow  for  the  lightest  possible 
configuration.  This  frees  up  ship’s  weight  restrictions  and  allows  for  more 
missiles  to  be  carried. 

•  Missiles  fire  vertically,  horizontally,  or  at  any  angle  in  between.  Thus, 
allowing  CCL’s  to  be  outfitted  in  geometrically  constraining  spaces  such  as 
the  narrow  hull  designs  of  the  swath  and  surface-effects  ships. 

•  A  shock  collar  mitigates  the  underwater  shock  at  the  main  deck  rather  than 
the  inner  bottom  or  platform  level. 

•  A  wide  assortment  of  different  weapons  loadouts,  such  as  ships  self 
defense  weapons  (Seasparrow  AAW  missiles),  electronic  warfare  devices 
(chaff  and  flares),  torpedoes  (MK  46  and  50),  remote  sensing  devices 
(sonobuoys)  in  addition  to  the  current  capabilities  of  the  MK  41  VLS 
(TLAM,  SLAM,  NATACMS,  SM-2,  etc.). 
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Figure  3.  Shock  Collar. 
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The  CCL,  although  different  from  the  MK  41  VLS,  is  not  intended  to  replace  it,  but 
merely  provides  ship  builders  with  a  new  set  of  options  for  launching  a  wide  variety  of 
weapons. 

A  CCL  module  sized  to  fit  in  the  after  VLS  space  of  a  DDG-5 1  ship  contains  up  to 
80  Tomahawk  missiles.  This  weapons  loadout  contains  approximately  24,000  lb.  of  solid 
propellant,  30,000  lb.  of  explosives,  and  35,300  lb.  of  liquid  fuel  [Ref.  1],  Thus  the  fires 
experienced  by  British  ships  in  the  Falklands  conflict  and  the  USS  STARK  (FFG  31) 
involving  this  much  high  energy  material  has  provided  a  vital  insight  into  the  dangers  and 
difficulties  of  fire  engulfing  the  spaces  surrounding  a  CCL  magazine.  These  incidents  have 
aided  greatly  in  the  development  of  realistic  scenarios  to  predict  ship  survivability. 

C.  PREVIOUS  WORK 

Bowman  and  Lee  [Ref.  1  ]  developed  a  one  dimensional  heat  transfer  model  to 
predict  the  time-temperature  profiles  in  the  missile  canisters  of  the  MK  41  VLS  system  for 
DLG-5 1  FLT  I,  II  and  IIA  ships.  Additionally,  this  model  was  used  to  evaluate  the 
effectiveness  of  high  temperature  thermal  insulation  to  reduce  the  heat  transfer  rates. 

Their  scenarios  used  a  fire  in  an  adjacent  compartment,  which  increased  temperatures 
from  ambient  conditions  to  2000°F  in  a  period  of  five  minutes,  and  then  maintained  that 
temperature  for  the  duration  of  the  fire. 

Mansfield  [Ref.  3]  conducted  an  analysis  of  ordnance  heating  in  steel  walled 
compartments.  He  developed  a  computer  algorithm  that  introduces  heat  into  the  fire 
compartment  according  to  any  selected  fire  size,  and  follows  the  heat  transfer  into  the 
magazine  and  other  compartments  that  surround  the  fire  compartment.  The  output  of  this 
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algorithm  is  the  calculation  of  the  temperature  of  the  fluid  medium  in  the  compartment  and 
the  surface  temperature  of  the  surrounding  bulkheads,  deck  and  overhead.  It  does  not, 
however,  take  into  account  any  varying  geometry  within  the  magazine. 

D.  OBJECTIVES 

The  purpose  of  this  study  is  to  model  a  particular  CCL  within  the  module, 
specifically  the  CCL  located  in  the  first  row  along  the  centerline  and  corner,  for  a  DDG-5 1 
FLT-IIA  ship.  This  one  dimensional  lumped  parameter  model  describes  the  time- 
temperature  profile  in  a  Tomahawk  missile  to  predict  the  time  it  takes  to  reach  a  critical 
value  or  “cook-off”  temperature  of  propellant  at  the  geometric  center  of  the  missile. 
Different  fire  intensities,  temperatures  and  heat  fluxes  will  be  used  to  provide  a  wide 
baseline  for  comparison  to  real  world  scenarios. 


6 


II.  MODEL 


A.  FORMULATION 

The  computer  algorithm  modeling  the  transfer  of  heat  from  a  compartment  on  fire 
to  a  missile  contained  in  a  Concentric  Canister  Launcher  (CCL)  magazine  is  based  on  the 
electrical  circuit  analogy  for  heat  transfer.  Just  as  an  electrical  circuit  is  associated  with 
the  conduction  of  electricity,  a  thermal  circuit  may  be  associated  with  the  transfer  of  heat 
[Ref.  4],  Therefore,  the  first  step  in  the  creation  of  the  model  is  to  develop  a  thermal 
network  describing  the  heat  transfer  from  the  space  to  the  missile  taking  into  account  the 
geometry  of  the  CCL  magazine.  Time  dependence  is  accounted  for  in  the  network  by 
assigning  a  capacitance  for  the  solid  material  of  the  bulkhead,  canister  walls,  and  missile. 
The  input  parameter  to  the  algorithm  is  the  rate  of  temperature  increase  on  the  outer 
surface  of  the  bulkhead  from  ambient  conditions  to  some  specified  steady  value  based  on 
the  fire  scenario  analyzed.  Thermal  conduction,  convection  and  surface  radiation  are  all 
accounted  for  in  the  model.  Each  one  will  be  discussed  along  with  the  associated 
correlations,  thermophysical  properties,  and  geometry  associated  with  its  position  in  the 
thermal  resistance  network.  The  computer  code  for  the  heat  transfer  model  is  written  in 
MATLAB  4.2.c.  A  schematic  of  the  thermal  resistance  network  is  given  in  Figure  4. 

1.  Conduction 

Conduction  takes  place  in  four  locations  in  the  resistance  network.  The  first  is  the 
bulkhead  which  is  modeled  as  a  plane  wall  constructed  of  mild  steel,  AISI  1020.  The 
thermophysical  property  data  for  this  material  were  obtained  from  Touloukian’s, 
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“Thermophysical  Properties  of  Matter”  [Ref.  6],  The  data  were  curve  fitted  using 
MATLAB  to  form  a  polynomial  to  develop  a  smooth  curve  that  “best  fits”  the  property 
data  and  can  be  evaluated  easily  by  the  algorithm.  The  properties  required  for  the 
algorithm  are  thermal  conductivity  (ksl),  specific  heat  (cst),  and  density  (pst).  Each  were 
evaluated  at  the  average  temperature  between  the  two  surfaces  of  the  bulkhead.  The 
thermal  resistance  for  conduction  through  the  bulkhead  is  given  by  the  equation 


cond 


Ax 
KtA  i 


(1) 


where  (Ax)  is  the  bulkhead  thickness  and  (A,)  is  the  area  normal  to  heat  transfer. 
Therefore,  the  thermal  energy  transferred  through  the  solid  due  to  conduction  is  described 
by  the  relationship; 
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The  capacitance  assigned  to  the  bulkhead  to  account  for  time  dependence  is  expressed  in 
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Figure  4.  Thermal  Resistance  Network  of  the  CCL  Magazine. 
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where  the  average  temperature  is  used  to  describe  the  temperature  profile  through  the 
solid.  This  model  assumes  temperature  profiles  to  be  linear  at  any  instant  through  solid 
material.  The  capacitance  equation  for  stored  energy  remains  the  same  for  the  cylinder 
walls  and  missile,  except  in  the  thermophysical  properties  associated  with  the  different 
material  used  in  their  construction  and  in  the  calculation  of  the  volume,  which  takes  into 
account  their  different  geometric  shapes. 

Conduction  through  the  outer  and  inner  cylinder  walls  is  evaluated  similar  to  the 
bulkhead,  except  that  it  is  modeled  in  cylindrical  coordinates  and  through  walls 
constructed  of  titanium.  Their  thermophysical  properties  are  again  obtained  from 
Touloukian,  and  are  curve  fitted  with  MATLAB  to  form  a  polynomial  for  use  by  the 
algorithm  [Ref.  6],  The  thermal  resistance  equations  for  conduction  in  the  outer  and  inner 
cylinder  wall  are; 
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The  missile’s  thermal  resistance  is  constructed  in  the  same  manner  as  above  except  the 
missile  thermophysical  properties  are  assumed  to  be  aluminum  and  uniform  throughout. 

2.  Convection 

Natural  convection  occurs  in  three  areas  in  the  resistance  network.  The  first  is  the 
space  between  the  bulkhead  and  the  outer  surface  of  the  CCL.  The  second  is  the  annulus 
between  the  inner  and  outer  cylinders,  which  function  as  the  gas  management  system 
during  missile  launch.  The  third  is  the  space  between  the  inner  cylinder  and  the  surface  of 
the  missile.  The  fluid  medium  of  each  area  is  assumed  to  be  air  and  is  evaluated  at  a  mean 
boundary  layer  temperature,  termed  the  film  temperature.  The  properties  required  for 
computing  the  convective  resistance  are  thermal  conductivity  (k^),  kinematic  viscosity 
(v),  volumetric  thermal  expansion  coefficient  (P),  and  Prandtl  number  (Pr).  The  thermal 
resistance  for  convection  for  each  air  space  is  given  by  the  equation 


R  = 


conv 


l 

h*A 


(6) 


where  the  calculation  of  the  average  convection  heat  transfer  coefficient  (  h  )  is  based 
upon  the  Nusselt  number  (  Nu  ),  thermophysical  properties  of  air,  and  either  the  height 
(L)  of  the  CCL  or  the  bulkhead. 

_  Wxk  . 

h=_u_mL 
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The  thermophysical  properties  of  air  were  curve  fitted  from  tables  in  Reference  4  for  use 
by  the  computer  algorithm.  The  Nusselt  number  is  a  dimension  less  heat  transfer 
coefficient  and  provides  a  measure  of  the  convection  heat  transfer  occurring  at  the  surface. 
Empirical  correlations  for  the  Nusselt  number  have  been  developed  based  on  the  geometry 
of  each  enclosed  air  space. 

The  first  air  space  Nusselt  numbers  are  calculated  with  respect  to  the  bulkhead  and 
the  CCL’s  free  convection  boundary  layer.  The  bulkhead’s  Nusselt  number  is  determined 
by  using  a  correlation  for  a  vertical  plate  developed  by  Churchill  and  Chu  that  may  be 
applied  over  an  entire  range  of  Rayleigh  numbers  [Ref.  4], 
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0.387/to/ 


lJ  °-492l 

9  ' 
16 

l  Pr  J 

(8) 


The  same  equation  may  be  applied  to  the  outer  surface  of  the  CCL  if  the  boundary  layer 
thickness  is  much  less  than  the  diameter  of  the  cylinder.  To  check  for  this  condition  the 
following  criterion  must  be  met  [Ref.  4], 
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can 
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If  this  inequality  is  not  met,  then  curvature  effects  are  appreciable  and  must  be  accounted 
for  using  the  following  correlation  by  Le  Fevre  and  Ede  [Ref  5]; 


7GrLPr2  U  4(272+315 Pr)L 
5(20  +21  Pr)  J  3  5(64 +63 Pr)d 


(10) 


The  Grashof  (Gr)  and  Rayleigh  (Ra)  numbers  used  throughout  these  correlations  are  given 
by: 


(ii) 


Ra^Grfr 


(12) 


After  each  Nusselt  correlation  is  computed,  the  average  convective  heat  transfer 
coefficient  is  calculated  from  Equation  7.  Assuming  both  surfaces  act  as  external  free 
convection  flows,  resistances  are  calculated  with  respect  to  the  boundary  layers  along  the 
bulkhead  and  CCL’s  surface  using  Equation  6.  These  resistances  are  added  in  series  to 
form  the  total  convective  resistance  between  the  bulkhead  and  the  outer  surface  of  the 
CCL.  A  schematic  of  the  resistances  is  given  in  Figure  5. 


13 


Figure  5.  1st  Air  Space  Resistance  Network. 


The  final  two  air  spaces  are  circular  annuli,  and  the  Nusselt  correlation  used  was 
developed  by  Churchill  and  Ozoe  [Ref.  7], 
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In  the  two  circular  annuli  only  one  resistance  is  calculated  to  describe  the 
resistance  across  the  enclosed  space.  The  resistance  and  average  convective  heat  transfer 
coefficient  is  again  calculated  using  Equation  6  and  7,  respectively.  A  schematic  of  the 
resistances  across  the  annuli  is  given  in  Figure  6. 


Figure  6.  Resistance  Network  Across  Annuli. 


3.  Radiation 

Radiation  affects  heat  transfer  in  the  same  air  spaces  as  convection.  Therefore,  the 
thermal  resistance  for  radiation  and  convection  act  in  parallel,  and  can  be  combined  to 
obtain  a  single  effective  surface  resistance  [Ref.4], 


15 
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The  total  resistance  to  radiation  exchange  between  the  surfaces  is  comprised  of  the  two 
surface  resistances  and  the  geometrical  resistance. 


Rr 


tot 


1 

£2^2  ^23^2  C3^3 


(16) 


The  radiative  property  used  in  the  resistance  calculations  is  emissivity  (e).  The  emissivity 
of  the  different  materials  is  obtained  from  Touloukian  [Ref.  6],  and  curve  fitted  to  form  a 
polynomial  for  use  in  the  computer  algorithm.  The  view  factor  (F23)  is  explained  later. 
Therefore,  the  net  radiation  exchange  between  surfaces  is  expressed  as; 


o(7-,4-r,4) 


Rr 


tot 


(17) 


However,  this  equation  is  non-linear  to  the  fourth  power,  and  the  overall  radiative 
resistance  (R^)  must  be  solved,  so  the  heat  rate  due  to  radiation  takes  the  form; 


Qrad 


V2-tj 

Rrad 


(18) 
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and  can  be  used  in  conjunction  with  the  convective  resistance.  This  is  accomplished  using 
two  different  methods.  The  first  is  to  linearize  the  equation  by  expanding  in  a  Taylor 
series,  and  the  second  is  to  solve  it  explicitly.  Using  a  Taylor  series  expansion  on 
Equation  17  takes  the  form; 


Qrad=(lo+' 


dq 

*r)T. 


(T-T,)+0[(T-T,f] 


(19) 


and  after  linearizing,  the  resistance  due  to  radiation  becomes 

1  _  90  toTl 

Rrad~(T2-T3)+  Rrtot 


(20) 


To  solve  for  the  overall  radiative  resistance  explicitly.  Equations  17  and  18  are  set  equal  to 
each  other  and  solved  algebraically.  The  overall  resistance  due  to  radiation  is  then 
expressed  as: 


1 

Rrad 


YT2*T^T2*T3)\ 


(21) 


Note  that  the  resistance  is  temperature  dependent.  These  two  processes  are  used  for  all 
radiative  resistances  in  the  thermal  resistance  network,  and  will  be  compared  in  the  final 
results.  The  difference  in  the  geometric  shapes  of  the  enclosures  are  taken  into  account  by 
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calculating  the  fraction  of  the  radiation  leaving  one  surface  and  intercepted  by  another. 
This  is  termed  the  view  factor  (F23).  In  the  first  air  space,  the  view  factor  is  calculated  as 
an  infinite  plane  to  a  row  of  cylinders  [Ref.  9], 


2  / 
F23=l  -D2)2 +D  arctan 
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2 

2 


(22) 


where 


(23) 


In  the  final  two  air  spaces,  the  view  factors  are  calculated  as  the  interior  of  an  outer  right 
circular  cylinder  of  finite  length  to  the  exterior  of  a  inner  right  circular  coaxial  cylinder  of 
finite  length  [Ref.  10], 
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where 
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(25) 


H= 


R= 


(26) 


The  thermal  energy  transferred  through  the  air  space  is  described  by  the  relationship 


q= 


(27) 


where  (R,)  is  the  total  resistance  due  to  radiation  and  convection,  Equation  15. 

The  radiative  resistance  network  thus  far  models  a  CCL  on  the  first  row  near  the 
centerline  of  the  ship.  As  mentioned  in  the  objectives,  the  CCL  in  the  comer  of  the 
magazine  is  also  modeled,  accounting  for  the  radiation  effects  of  the  adjoining  bulkhead. 
The  adjoining  bulkhead  is  modeled  as  a  reradiating  surface.  Equation  16  is  modified  to 
account  for  the  radiation  leaving  the  effected  bulkhead,  reflected  off  the  adjoining 
bulkhead  and  absorbed  by  the  surface  of  the  CCL,  and  takes  the  form; 
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The  view  factor  from  the  reradiating  bulkhead  to  the  canister  (FR3)  is  the  same  as 
the  view  factor  from  the  effected  bulkhead  to  the  canister  (F23),  Equation  22.  However,  a 
new  view  factor  (F2R)  to  determine  the  amount  of  radiation  received  by  the  adjoining 
bulkhead  is  calculated  as  perpendicular  plates  with  a  common  edge  [Ref.  4], 


F  = 

1  1R 


(w2) 

i+ 

(w2y 

(!) 

— 

— 

— 

KJ 

V  WRJ 

2 


(29) 


These  equations  are  solved  as  previously  discussed.  A  schematic  of  the  surface  radiation 


resistances  is  given  in  Figure  7. 


Figure  7.  1st  Air  Space  Resistance  Network  (Reradiating  Surface). 
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4.  Finite  Difference  Equations 

The  capacitance  and  transferred  thermal  energy  from  conduction  and  the  effects  of 
convection  and  radiation  are  combined  to  form  a  set  of  explicit  finite  difference  equations 
to  solve  for  the  temperature  at  each  node  in  the  thermal  resistance  network.  Each  is 
solved  simultaneously  at  each  instant  in  time,  which  is  defined  by  the  time  step  specified  in 
the  algorithm.  These  finite  difference  equations  were  derived  by  conducting  an  energy 
balance  about  each  temperature  node.  Because  this  is  a  time  dependent  model  due  to  the 
capacitance  in  the  solid  material  of  the  bulkhead,  canister  walls  and  missile,  the  thermal 
energy  stored  in  the  solid  material  is  taken  into  account  in  the  energy  balance. 

^in  stored  (30) 


All  heat  flow  is  assumed  to  flow  into  the  node,  and  that  each  node  in  the  resistance 
network  is  initially  at  ambient  conditions  (300  K).  The  result  of  the  energy  balance  yields 
equations  for  the  nodal  temperature  at  a  time  increment. 
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T7(p+ 1)=2 
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B.  LIMITATIONS 

In  the  construction  of  this  model,  a  high  level  of  accuracy  was  not  anticipated  due 
to  its  one  dimensionality.  However,  the  goal  is  to  develop  a  model  that  achieves 
moderately  accurate  solutions  to  provide  baseline  data  on  the  cook-off  of  Tomahawk 
missiles  contained  in  a  CCL  magazine.  The  primary  limitations  and  associated  errors  with 
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this  model  are  in  the  finite-differencing  method,  thermophysical  property  data,  geometric 
assumptions,  and  simplification  of  the  heat  transfer  process. 

1.  Finite  Difference  Method 

The  finite  difference  solution  used  is  the  explicit  method.  In  this  solution  method, 
the  unknown  nodal  temperatures  for  the  new  (p+1)  time  are  determined  from  known 
nodal  temperatures  at  the  previous  (p)  time.  The  accuracy  of  this  method  depends  on  the 
distance  between  nodes  (Ax)  and  the  time  step  (At),  as  these  values  decrease  the  accuracy 
of  the  solution  increases.  The  distance  between  nodes  is  chosen  based  on  the  geometry  of 
the  model,  and  the  time  step  is  chosen  based  on  the  stability  requirements.  The  difficulty 
with  the  explicit  method  is  that  it  is  not  unconditionally  stable.  This  condition  is 
characterized  by  oscillations  which  diverge  from  a  steady  value  instead  of  converging  as 
anticipated  in  transient  conduction  problems.  To  determine  the  time  step  at  which  the 
algorithm  used  in  the  model  is  stable,  the  following  criterion  for  a  one  dimensional  node  is 
solved  for  At. 


a  At  1 
(Ax)2  2 


(38) 


The  limiting  time  step  in  this  model  is  governed  by  the  inner  cylinder  wall  based  on  its 
values  of  wall  thickness  (Ax)  and  thermal  diffiisivity  (a).  Therefore,  maximum  allowable 
time  step  is  0.37  seconds.  The  more  the  time  step  is  reduced  below  0.37  seconds,  the 
more  accurate  the  solution.  However,  the  number  of  iterations  to  achieve  the  desired  time 
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span  increases  rapidly  and  makes  it  computationally  intensive.  Therefore,  0.37  seconds  is 
used  in  all  calculations  to  optimize  the  computing  time. 

2.  Property  Data 

The  thermophysical  properties  obtained  from  Touloukian  [Ref.  6],  which  are 
entering  arguments  throughout  the  algorithm,  are  a  large  source  of  error.  Not  only  is 
there  a  error  associated  with  curve  fitting  of  the  data,  but  also  in  the  experimental 
procedure  in  obtaining  the  properties.  The  error  reported  by  Touloukian  [Ref.  6]  in  the 
experimental  procedures  used  to  obtain  the  thermophysical  property  data  of  the  substances 
used  in  the  CCL  are  listed  in  Table  1. 


EXPERIMENTAL  ERROR 

STEEL 

TITANIUM 

AIR 

Specific  Heat  (cp) 

4% 

0.2% 

n/a 

Thermal  Conductivity  (k) 

1% 

5-15% 

1-5% 

Emissivity  (e) 

2% 

8% 

n/a 

Thermal  Diffusivity  (a) 

5-10% 

10% 

n/a 

Kinematic  Viscosity  (v) 

n/a 

n/a 

2% 

Table  1 .  Experimental  Error  in  the  Measurement  of  Thermophysical  Properties. 


In  order  for  the  algorithm  to  evaluate  the  properties  of  the  metals  and  air,  a  curve 
fit  or  regression  was  conducted  on  the  data.  In  this  procedure,  a  curve  was  developed  that 
‘best  fits”  the  data  but  does  not  necessarily  pass  through  the  data  points.  The  “best  fit”  or 
least  squares  curve  fitting  seeks  to  minimize  the  sum  of  the  squared  error.  The  error  is  the 
distance  between  the  polynomial  curve  and  the  actual  data  point.  Squaring  this  distance  at 
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each  data  point  and  adding  the  squared  distances  together  is  the  sum  of  the  squared  error. 
Therefore,  the  curve  generated  by  the  polynomial  is  the  curve  that  makes  the  sum  of  the 
errors  as  small  as  it  can  be,  “best  fit”  [Ref  8],  The  error  in  relation  to  the  data  and  the 
polynomial  generated  by  the  curve  fit  was  less  than  ten  percent  in  all  cases,  which  is  within 
the  range  of  the  experimental  data  error. 

3.  Geometry 

In  the  modeling  of  the  CCL  each  geometric  shape  is  considered  to  be  solid 
polished  material  with  no  attachments.  Obvious  architecture  with  respect  to  a  functioning 
weapon  system  have  been  neglected  such  as  wire  ways,  electronic  packaging,  supporting 
frame  work,  etc. 

Several  assumptions  have  been  made  with  respect  to  radiation.  The  metal  on  all 
surfaces  throughout  the  model  are  assumed  to  be  unpainted  and  polished.  Thus,  the 
emissivity  is  assumed  to  be  much  higher  than  that  of  a  painted  surface.  Additionally,  the 
model  accounts  only  for  the  heat  emitted  from  the  bulkhead  to  the  modeled  CCL  and 
neglects  radiation  effects  from  the  surrounding  canisters.  Finally,  radiation  with  respect  to 
the  other  walls  in  the  room  on  fire  to  the  CCL  bulkhead  is  neglected,  as  well  as,  the 
convective  resistance  on  the  outer  surface  of  the  bulkhead. 

The  final  assumption  in  the  model  is  that  the  missile  has  the  properties  of  aluminum 
and  is  solid  throughout.  The  missile  electronics  and  cabling  are  neglected.  Aluminum  was 
chosen  because  it  is  used  consistently  throughout  the  missile.  Its  property  data  would 
provide  a  good  assumption  as  to  how  the  missile  reacts  to  the  heat  energy  absorbed  and 
provide  a  temperature  profile  at  the  center  of  the  missile. 
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III.  FIRE  SCENARIOS 


To  develop  accurate  time-temperature  profiles  describing  the  heat-up  rate  of  the 
CCL,  credible  fire  scenarios  are  identified  based  on  the  type  of  combustible  material 
fueling  the  fire  and  its  probable  intensity.  Four  scenarios  are  analyzed  by  the  CCL  heat 
transfer  model  to  achieve  a  good  understanding  of  how  the  CCL’s  heat-up  rate  will  differ 
when  exposed  to  different  types  of  fires. 

The  first  scenario  analyzed  consisted  of  a  compartment  fire  involving  common 
shipboard  fuels  (JP-5,  JP-4,  F-76,  etc).  The  time-temperature  profile  is  characterized  by  a 
increase  in  temperature  from  ambient  conditions  to  2000°F  (1366  K)  over  a  period  of  five 
minutes  and  is  then  maintained  at  the  maximum  temperature  for  the  duration  of  the  fire 
[Ref.  1].  This  scenario  is  typically  used  to  evaluate  the  structural  integrity  of  the  effected 
space  both  during  and  after  the  fire  [Ref.  11], 

The  second  fire  scenario  is  typically  used  in  modeling  fires  fueled  by  common 
combustible  material,  such  as  a  berthing  compartment  fire.  The  time-temperature  profile 
is  characterized  by  a  increase  in  temperature  from  ambient  conditions  to  1700°F  (1200  K) 
over  a  period  of  fifteen  minutes,  and  then  maintains  a  steady  temperature  rise  of  84°F/hour 
(46K/hour)  [Ref.  11]. 

A  fire  fueled  by  a  missile’s  residual  solid  propellant,  as  experienced  by  the  USS 
STARK  (FFG  31),  is  the  most  intense  and  is  the  next  scenario  analyzed.  Propellant  fires 
are  normally  short  lived  but  can  reach  temperatures  in  excess  of  3000°F  (1922  K).  For  the 
purpose  of  this  scenario,  it  is  assumed  that  maximum  temperature  is  reached  in  thirty 
seconds,  and  that  the  temperature  is  then  maintained  until  the  self-ignition  temperature  of 
the  missile’s  solid  propellant  is  reached. 
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The  last  scenario  analyzed  was  developed  by  the  U.  S.  Navy’s  Insensitive 
Munitions  (IM)  Program  to  characterize  the  thermal  response  of  munitions  exposed  to 
slow  cook-off  conditions.  The  slow  cook-off  condition  is  described  by  a  temperature  rise 
of  6°F/hour  (3.33  K/hour)  [Ref.  1], 

The  IM  Program  also  developed  a  fire  scenario  to  describe  fast  cook-off,  that 
immerses  the  missile  in  a  JP-5  pool  fire  in  which  the  average  temperature  must  be  at  or 
above  1600°F  [Ref.  1],  This  scenario  is  most  likely  on  ordnance  stored  on  a  flight  deck 
rather  than  a  missile  contained  in  a  CCL  magazine.  Therefore,  the  fast  cook-off  scenario 
is  not  analyzed. 

The  test  data  from  the  Insensitive  Munitions  program  shows  that  the  missile’s  solid 
propellant  motor  (MK  106  Booster)  and  explosives  will  self-ignite  between  300  and  400°F 
(422  and  478  K),  and  the  liquid  propellant  motor,  fueled  by  JP-10  or  RJ-4,  can  self-ignite 
at  temperatures  as  low  as  460°F  (511  K).  Depending  on  the  rate  of  temperature  increase, 
self-ignition  of  the  energetic  material  can  range  from  a  slow  burn  to  an  immediate 
deflagration.  [Ref.  1] 

Table  2  summarizes  the  aforementioned  fire  scenarios  and  assigns  a  scenario 
number.  This  number  is  used  to  refer  to  the  various  fire  scenarios  throughout  the  rest  of 
the  text. 
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FIRE  SCENARIOS 


MAXIMUM 
TEMPERATURE  (K) 

RATE  OF  TEMP 
INCREASE 

1366 

3.55  K/second 

n/a 

80K/min<  1186K 

0.77  K/min  >  1 186  K 

1922 

54.07  K/second 

n/a 

3.33  K/hour 

Table  2.  Fire  Scenarios. 
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IV.  RESULTS 


As  mentioned  in  the  objectives,  this  study  models  two  CCLs  within  the  ship’s 
weapon  module.  The  first  is  located  in  the  first  row  along  the  centerline  of  the  ship,  and 
the  second  is  the  corner  CCL  in  which  the  adjacent  wall  is  considered  a  reradiating 
surface.  The  four  fire  scenarios  mentioned  in  the  previous  chapter  are  analyzed  for  each 
CCL  location.  Additionally,  each  scenario  is  conducted  twice.  The  first  solves  for  surface 
radiative  resistance  explicitly  and  the  second  with  a  Taylor  series  expansion.  After  the 
scenario  number,  the  explicit  method  and  Taylor  series  expansion  methods  will  be  denoted 
by  the  letter  (A)  and  (B),  respectively. 

The  temperature  range  in  which  the  missile’s  components,  specifically  the  warhead 
and  solid  propellant  motor,  will  self-ignite  is  between  422  and  478  K.  The  liquid 
propellant  will  self-ignite  as  low  as  5 1 1  K.  The  times  to  reach  these  temperatures  for  the 
center  and  comer  CCL  are  listed  in  Table  3  and  4,  respectively.  Time-temperature  profiles 
are  plotted  in  Figures  8  through  15  for  the  center  CCL  and  Figures  16  through  23  for  the 
comer  CCL. 

The  two  methods  used  to  determine  the  surface  radiation,  described  in  Chapter  II, 
are  compared  because  of  the  large  differences  experienced  in  each  scenario’s  time- 
temperature  profiles.  The  Taylor  series  expansion  produces  two  types  of  error,  truncation 
and  round-off.  The  truncation  error  is  due  to  the  approximation  formula.  In  general, 
increasing  the  order  of  the  Taylor  series  reduces  the  truncation  error  but  increases  the 
round-off  error.  For  this  model  both  types  of  errors  are  probably  high  for  two  reasons: 
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the  expansion  was  only  carried  out  to  the  first  derivative  to  reduce  the  amount  of 
computations,  and  the  round  off  error  is  high  due  to  the  curve  fitted  thermophysical 
properties  used  throughout  the  model.  Solving  for  the  surface  radiative  resistance 
explicitly  eliminates  the  truncation  error,  however,  it  is  still  affected  by  the  round  off  error. 
Solving  for  the  surface  radiative  resistances  with  a  Taylor  series  expansion  produces  a 
more  conservative  time-temperature  profile,  and  is  the  recommended  method. 
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TIME  TO  COOK-OFF 

SCENARIO  # 

TIME  TO  REACH 
511  K  (min) 

TIME  TO  REACH 
422-478  K  (min) 

ASSOCIATED 

FIGURE 

1A 

204.42 

141.06-181.02 

8 

IB 

152.26 

108.48-136.32 

9 

2A 

254.40 

187.52-230.96 

10 

2B 

198.68 

147.43-180.73 

11 

3A 

64.85 

49.71-59.28 

12 

3B 

53.25 

42.23-49.25 

13 

4A 

inf 

inf 

14 

4B 

inf 

inf 

15 

Table  3.  Cook-Off  Times  for  Scenarios  (Center  CCL). 


TIME  TO  COOK-OFF 

SCENARIO  # 

TIME  TO  REACH 
511  K  (min) 

TIME  TO  REACH 
422-478  K  (min) 

ASSOCIATED 

FIGURE 

1A 

199.88 

137.76-176.95 

16 

IB 

149.59 

106.43-133.88 

17 

2A 

250.04 

183.92-226.88 

18 

2B 

195.83 

145.02-178.03 

19 

3A 

63.50 

48.67-58.06 

20 

3B 

52.42 

41.56-48.48 

21 

4A 

inf 

inf 

22 

4B 

inf 

inf 

23 

Table  4.  Cook-Off  Times  for  Scenarios  (Comer  CCL). 
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Figure  8.  Scenario  lA-Missile  Heat-Up  Rate. 
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Figure  9.  Scenario  IB-Missile  Heat-Up  Rate. 
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Figure  10.  Scenario  2A-Missile  Heat-Up  Rate. 
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Figure  13.  Scenario  3B-Missile  Heat-Up  Rate. 
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Figure  16.  Scenario  1  A-Missile  Heat-Up  Rate  w/  Reradiating  Surface. 
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Figure  20.  Scenario  3  A-Missile  Heat-Up  Rate  w/  Reradiating  Surface. 
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Figure  21 .  Scenario  3B-Missile  Heat-Up  Rate  w/  Reradiating  Surface. 


47 


Temperature  (K) 


8 


Temperature  (K) 


Figure  23.  Scenario  4B-Missile  Heat-Up  Rate  w/  Reradiating  Surface. 
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V.  CONCLUSIONS 


The  fire  scenarios  analyzed  by  the  CCL  heat  transfer  model  in  this  study  showed 
that  the  CCL  can  withstand  an  intense  conflagration  if  the  bulkheads  surrounding  the 
magazine  remain  intact.  In  the  most  intense  fire  situation,  Scenario  3,  ship’s  personnel 
have  approximately  forty-five  minutes  to  gain  control  of  the  fire  or  cool  the  weapons  in 
the  magazine  prior  to  cook-off.  Conceivably,  automated  fire  protection  equipment  could 
increase  the  amount  of  time  firefighters  have  to  control  the  fire  without  the  immediate 
hazard  of  cook-off 

In  its  current  stage  of  development,  the  CCL’s  installed  fire  protection  features 
have  not  been  designed.  However,  the  uniqueness  of  the  CCL  module  will  require  some 
changes,  as  compared  to  the  fire  protection  equipment  currently  installed  in  the  MK  41 
VLS.  The  MK  41  VLS  fire  protection  system  consist  of  a  deluge  system  installed  in  the 
canister  designed  to  cool  the  missile  in  the  event  of  an  inadvertent  motor  ignition  or 
extreme  over  temperature  situation,  and  a  sprinkling  system  designed  to  provide  cooling 
external  to  its  canisters  should  it  sense  a  rapid  temperature  rise  [Ref.  1],  Because  the  CCL 
acts  both  as  a  launcher  and  transport  device,  introducing  a  fire  protection  system  into  the 
CCL  will  hamper  its  primary  advantage  of  being  a  totally  integrated  package. 

A  possible  alternative  to  the  current  deluge  system  would  be  a  similar  system  of 
two  or  more  mounted  spray  nozzles  directed  to  impinge  the  canister  at  the  top  outer 
surface  of  the  CCL  providing  a  cooling  film  around  the  entire  canister  but  not  wetting  the 
missile.  In  addition  to  the  outer  deluge  system,  a  magazine  sprinkler  should  also  be 
installed  to  cool  the  entire  magazine  should  temperatures  increase  too  rapidly.  The 
advantages  of  keeping  the  fire  protection  system  outside  of  the  CCL  are  two  fold; 
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•  If  there  is  an  inadvertent  sprinkling  of  the  magazine  or  activation  of  the 
deluge  system,  the  missile  is  not  effected  by  the  water  and  the  costly 
process  of  sending  the  weapon  to  a  rework  facility  for  repairs  and 
recertification  or  decertification  and  dismantling  is  avoided.  It  is  possible 
that  after  a  incident  of  this  type  onboard  data  transmission,  continuity  and 
accuracy  checks  can  be  conducted  to  verify  the  missile  is  fully  functional 
after  the  module  and  electrical  connections  are  cleaned. 

•  The  second  advantage,  and  most  important,  is  to  maintain  the  ease  of  the 
installation  and  removal  process  from  the  ship’s  weapon  module.  Installing 
piping  connections  for  a  deluge  system,  which  of  course  will  have  to  meet 
high  shock  criteria,  will  make  a  relatively  easy  installation/removal 
operation  much  more  cumbersome. 

The  disadvantage  of  water  systems  is  the  generation  of  steam  in  intense  fires 
resulting  in  a  steam  explosion.  There  are  other  types  of  systems,  such  as  inert  gas  fire 
protection  systems,  which  may  be  worthy  of  investigation.  But  when  dealing  with  intense 
fires  with  long  durations,  as  in  the  USS  STARK,  the  availability  and  cooling  effects  of 
seawater  proved  to  be  effective.  Investigating  the  installation  of  a  relief  system  should 
steam  accumulate  in  the  module  during  a  fire  is  recommended. 

Another  characteristic  of  missile  solid  propellent,  brought  to  our  attention  by 
NAVAIR,  China  Lake,  which  warrants  mentioning,  is  that  prior  to  reaching  self-ignition 
temperature,  the  propellant  begins  to  off-gas  some  of  its  volatile  materials  [Ref.  11],  The 
temperature  at  which  this  occurs  is  yet  to  be  identified.  From  the  standpoint  of  any 
exposure  to  an  over-temperature  situation,  a  missile  may  be  damaged  even  if  cook-off  did 
not  occur.  However,  an  interesting  side  effect  was  noticed  that  if  the  solid  propellant  was 
remixed  after  the  off  gassing,  it  burned  longer  and  with  less  thermal  energy.  Follow  on 
research  is  recommended  to  determine  at  what  temperature  the  solid  propellent  begins  to 
off-gas. 
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In  general,  this  study  will  provide  CCL  project  managers  and  designers  with 
baseline  time-temperature  profiles  to  make  decisions  concerning  the  construction  and 
survivability  of  the  CCL  launching  system,  especially  with  respect  to  installed  fire 
protection  systems.  Additionally,  it  will  provide  ship’s  personnel  with  a  conservative 
missile  cook-off  time  based  on  the  type  of  fire,  so  that  the  damage  control  situation  can  be 
assessed  and  priorities  set. 
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VI.  RECOMMENDATIONS 


In  any  continuation  of  this  study  the  following  are  recommended: 

•  The  CCLs  modeled  are  sized  for  Tomahawk  missiles  and  are  the  same 
throughout  the  ship’s  weapon  module.  Future  studies  should  model  the 
ship’s  weapon  module  for  an  actual  deployment  weapon  configuration  in 
which  different  types  of  CCLs  carry  a  variety  of  weapons  in  support  of  all 
the  ship’s  mission  areas. 

•  Build  a  user  interface  at  the  MATLAB  command  window  to  easily  vary  the 
fire  scenario.  Currently,  to  vary  the  fire  scenarios  the  computer  code  must 
be  modified. 

•  Account  for  surface  radiation  from  other  CCLs  in  the  magazine. 

•  Utilize  temperature  data  from  fire  scenarios  currently  being  conducted  in 
the  ex-US  S  SHAD  WELL.  This  should  provide  more  accurate  data  and  a 
larger  variety  of  fire  curves  to  be  evaluated  by  the  CCL  Fire  Model. 
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APPENDIX  A.  CCL-FIRE  PROGRAM  CODE  (2A-CENTER) 

This  MATLAB  code  solves  for  radiation  explicitly  as  described  in  the  radiation 


section  of  Chapter  II.  This  code  is  set  up  for  the  center  CCL  and  scenario  2. 


format  bank 
clear 

%  DIMENSIONS  FOR  CANISTERS  &  SURROUNDING  VLS  WALL. 


w=l; 

% 

Lw=l; 

% 

Lc=l; 

% 

Lm=l; 

% 

delx=0.0048; 

% 

r=0. 00159; 

% 

ro=0.2586; 

% 

rl=0.2635; 

% 

r2=0.2667; 

% 

r3=0.3492; 

% 

r4=0.3571; 

% 

s=0.8382; 

% 

g=9.81; 

% 

sig=5.67e-008; 

% 

Al=w*Lw; 

% 

A2=A1; 

A3=2*pi*r4*Lc; 

% 

A4=2*pi*r3*Lc; 

% 

A5=2*pi*r2*Lc; 

% 

A6=2*pi*rl*Lc; 

% 

A7=2*pi*ro*Lc; 

% 

Vl=delx*w*Lw; 

V2=pi*Lc*(r4A2-r3A2); 

V3=pi*Lc*(r2A2-rlA2); 

V4=pi*roA2*Lm; 

[m]  Width  of  for  and  aft  bulkhead. 

[m]  Height  of  wall  surrounding  CCL  module. 

[m]  Height  of  inner/outer  canister  wall. 

[m]  Length  of  missile. 

[m]  Thickness  of  wall  surrounding  CCL  module, 
[m]  Near  center  of  missile  (1/16  in). 

[m]  Radius  of  missile. 

[m]  Inside  radius  of  inner  cylinder  wall. 

[m]  Outside  radius  of  inner  cylinder  wall. 

[m]  Inside  radius  of  outer  cylinder  wall. 

[m]  Outside  radius  of  outer  cylinder  wall. 

[m]  Distance  between  cylinder  CL. 

[m/sA2]  Gravitational  constant. 

[W/mA2*KA4]  Stefan-Boltzmann  constant. 
[mA2]Area  of  VLS  wall  perp  to  dir  of  heat  flux. 

[mA2]  Outer  surface  area  of  1st  cylinder. 

[mA2]  Inner  surface  area  of  1st  cylinder. 

[mA2]  Outer  surface  area  of  2nd  cylinder. 

[mA2]  Inner  surface  area  of  2nd  cylinder. 

[mA2]  Surface  area  of  missile. 


%  [K]  INITIAL  AMBIENT  CONDITIONS  TO  START  ITERATION. 

T2(l)=300; 

T3(l)=300; 
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T4(l)=300; 

T5(l)=300; 

T6(l)=300; 

T7(l)=300; 

T8(l)=300; 


%  DEFINES  TIME  STEP  AND  NUMBER  OF  ITERATIONS 

S=50000;  %  Number  of  iterations. 

delt=0.36;  %  [s]  Time  step. 


%  BEGINS  LOOP  TO  SOLVE  FOR  ALL  TEMPS  IN  RESISTANCE  NETWORK. 


for  p=l:S 
pl=p 


%  [s]  Time. 
t(p)=(delt*(p-l»; 

dT(p)=t(p)*  1.3333; 
dTA(p)=t(p)*0.0128; 

T 1  (p)=3 00+dT  (p); 

T 1  A(p)= 1 200+dT  A(p); 

if  Tl(p)<=1200 

Tl(p)=Tl(p); 

else 

Tl(p)=TlA(p); 

end 

t(p+l)=(delt*(p)); 
dT(p+l)=t(p+l)*  1.3333; 
dTA(p+l)=t(p+l)*0.0128; 

T 1  (p+ 1 )=3 00+dT  (p+ 1 ); 

T 1  A(p+ 1 )= 1 200+dT  A(p+ 1 ); 

if  Tl(p+l)<=1200 
Tl(p+l)=Tl(p+l); 
else 

Tl(p+l)=TlA(p+l); 

end 
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%%%  SOLVES  SECTION  1  [T2]  %%% 


%  [K]  Film  temperature  to  evaluate  properties  of  mild  steel 
Tfl(pHTl(p)+T2(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  bulkhead. 

kw  1  (p)=abs((  1 ,8869e-005*Tfl  (p)A2)-(6.9896e-002*Tfl  (p))+83 .43  8); 

%  [J/kg*K]  Specific  Heat  of  bulkhead. 

cpl(p)=abs(((-2.2828e-006)*Tfl(p)A3)+((4.6812e-003)*Tfl(p)A2)-(2.3469*Tfl(p))+800. 

99); 

%  [mA2/s]  Thermal  diffusivity  of  bulkhead. 

alfl  (p)=ab  s(((4 .  0 1 1 8e-0 1 4)  *  Tfl  (p)A3  )-((7 .  93  3  7e-0 1 1 )  *Tfl  (p)  A2)+(3 ,3733e-008*Tfl  (p)) 
+1.0388e-005); 

%  [kg/mA3]  Density  of  bulkhead, 
rho  1  (p)=kwl(p)/(cp  1  (p)*alfl  (p)); 

%  Solves  for  the  conductive  resistance  in  bulkhead. 

Rcl(p)=delx/(kwl(p)*A2); 

%  [K]  Film  temperature  to  evaluate  properties  in  1st  air  space. 

Tf2(p)=(T2(p)+T3(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  air. 

ka2(p)=abs((4.9198e-011*Tf2(p)A3)-(1.608e-007*Tf2(p)A2)+(2.014e-004*Tf2(p))-0.020 

796); 

%  [mA2/s]  Kinematic  viscosity  of  air. 

nu2(p)=abs(((7.786e-005*Tf2(p)A2)+(4.4053e-002*T£2(p))-2.4972)*1.0e-006); 

%  Prandtl  #. 

Pr2(p)=abs((-2.5065e-01 1  *Tf2(p)A3)+(8.4944e-008*Tf2(p)A2)-(l  .0417e-004*Tf2(p))+0. 
74553); 


%  [KA-1]  Exp.  Coefficient  of  air . 

B2(p)=l/Tf2(p); 

%  Emissivity  of  mild  steel. 

E2(p)=abs((-9 . 84 1 6e-008 *  T2(p)A2)+(2. 69 1 2e-004 *  T2(p))+0 .73699); 
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%  Emissivity  of  titanium. 

E3  (p)=abs((-5 .041 2e-008 *  T3  (p)  A2)+(2 . 4949e-004 *T3  (p))+0. 0 515 04); 

%  Solves  for  Rayleigh  and  Grashof  #'s  WRT  bulkhead  inner  surface. 
Gt2(p)=((g*B2(p)*(T2(p)-Tf2(p))*(LwA3))/(nu2(p)A2)); 

Ra2(p)=Gr2(p)*Pr2(p); 

%  Correlation  to  solve  for  the  average  Nusselt  #  over  a 
%  vertical  flat  plat.  Recommended  by  Churchill  and  Chu,  it  is 
%  good  over  entire  range  of  Rayleigh  #'s. 

Nu2(p)=(0.825+((0.387*(Ra2(p)A(l/6)))/((l+((0.492/Pr2(p))A(9/16)))A(8/27))))A2; 


%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  over  the  bulkhead  inner  surface. 
h2(p)=Nu2(p)*ka2(p)/Lw; 

Rco  1  a(p)=  1  /(h2(p)  *  A2); 


%  Solves  for  Rayleigh  and  Grashof  #'s  WRT  cylinder  outer  surface. 

Gr  3  (p)=(g*B2(p)*  (T  f2(p)-T3  (p))  *  (LcA3  ))/ (nu2(p)A2); 

Ra3  (p)=Gr3  (p)*Pr2(p); 

%  IF-THEN  LOOP  TO  DETERMINE  IF  THE  OUTER  WALL  OF  THE  CYLINDER 
%  CAN  BE  APPROXIMATED  AS  A  PLANE  WALL.  IF  A  PLANE  WALL 
%  CHURCHILL 

%  AND  CHU'S  CORRELATION  FOR  A  VERTICAL  SURFACE  IS  USED.  IF 
%  THERE  ARE  APPRECIABLE  CURVATURE  EFFECTS  A  CORRELATION 
%  DEVELOPED 

%  BY  LE  FEVRE  AND  EDE  BASED  ON  CYLINDER  HEIGHT. 


if  Gr3(p)==0 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)A2)/(5*(20+21*Pr2(p))))A0.25)+((4*Lc*(272+315*Pr2 

(p)))/(35*2*r4*(64+63  *Pr2(p)))); 

else 

if  ((2*r4)/Lc)>=3  5/(Gr3(p)Al/4) 

Nu3(p)=(0.825+((0.387*(Ra3(p)A(l/6)))/((l+((0.492/Pr2(p))A(9/16)))A(8/27))))A2; 

else 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)A2)/(5*(20+21*Pr2(p))))A0.25)+((4*Lc*(272+315*Pr2 

(P»)/(3  5  *  2  *  r4*  (64+63  *Pr2(p)))); 

end 

end 
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%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  over  outer  surface  of  the  cylinder. 
h3  (p)=(Nu3  (p)*ka2(p))/Lc; 

Rco  1  b(p)=  1  /(h3  (p)  *  A3  ) ; 

%  Adds  convective  resistances  in  series. 

Rco  1  (p)=Rco  1  a(p)+Rco  1  b(p); 

%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  infinite  plane  to  a  row  of  pipes  [C-6], 

D=(2*r4)/s; 

F23=l-((1  -D  A2)  A0 . 5)+(D*  atan(((  1  -DA2)/DA2)A0. 5)); 

%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

Rr  1  (P)=((  1  -E2(p))/ (E2(p)*  A2))+(  1  / ( A2  *F23))+((  1 -E3  (p))/ (E3  (p)*  A3  )); 

Rradl(p)=  1  /  ((sig/Rr  1  (p))*((T2(p)+T3  (p))  *  (T2(p)A2+T3  (p)  A2))); 

%  Solves  for  total  resistance,  [T3],  and  begins  next  iteration. 

Rt(p)=(Rco  1  (p)*Rrad  1  (p))/(Rrad  1  (p)+Rco  1  (p)); 

%  Finite  difference  equation— explicit  method. 

H 1  (p)=delt/ (rho  1  (p)  *  cp  1  (p)  *  V 1 ) ; 

T2(p+ 1  )=2  *  (((H 1  (p)/Rc  1  (p))  *(T  1  (p)-T2(p)))+((H  1  (p)/Rt(p))*  (T3  (p)-T2(p)))+((T2(p)+T 
l(p))/2))-Tl(p+l); 

if  T2(p+1)<300 

T2(p+1)=300; 

else 

T2(p+l)=T2(p+l); 

end 


%%%  SOLVES  SECTION  2  [T3]  %%% 

%  [K]  Film  temperature  to  evaluate  properties  of  titanium. 

Tf3(p)=(T3(p)+T4(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  canister  wall. 

kc3(p)=abs((  1.701  le-008*Tf3(p)A3)-(4.33 16e-005*Tf3(p)A2)+(2.7279e-002*Tf3(p))+20 
358); 
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%  [J/kg*K]  Specific  Heat  of  canister  wall. 

cp3(p)=abs(((4.6480e-007)*TD(p)A3)-((1.4813e-003)*Tf3(p)A2)+(1.5387*Tf3(p))+139.9 

3); 

%  [mA2/s]  Thermal  diffiisivity  of  canister  wall. 

alf3(p)=abs(((-4.532e-015)*Tf3(p)A3)+((1.7281e-011)*Tf3(p)A2)-(1.929e-008*TB(p))+l 

,3594e-005); 

%  [kg/mA3]  Density  of  canister  wall. 
rho3(p)=kc3(p)/(alf3(p)*cp3(p)); 

%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 

Rc2(p)=log(r4/r3  )/(2  *  pi  *  kc3  (p)  *  Lc) ; 

%  Finite  difference  equation-explicit  method. 

H2(p)=delt/(rho3  (p)  *cp3(p)*V2); 

T3(p+l)=2*(((H2(p)/Rt(p))*(T2(p)-T3(p)))+((H2(p)/Rc2(p))*(T4(p)-T3(p)))+((T3(p)+T 

4(p))/2))-T4(p); 


%%%  SOLVES  SECTION  3  [T4]  %%% 

%  [K]  Film  temperature  to  evaluate  air  properties  in  1st  annulus. 

Tf4(p)=(T4(p)+T5(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  air. 

ka4(p)=abs((4.9198e-011I|cTf4(p)A3)-(1.608e-007*Tf4(p)A2)+(2.014e-004*Tf4(p))-0.020 

796); 

%  [mA2/s]  Kinematic  viscosity  of  air. 

nu4(p)=abs(((7.786e-005*Tf4(p)A2)+(4.4053e-002*Tf4(p))-2.4972)*1.0e-006); 

%  Prandtl  #. 

Pr4(p)=abs((-2.5065e-011*Tf4(p)A3)+(8.4944e-008*Tf4(p)A2)-(1.0417e-004*Tf4(p))+0. 

74553); 

%  [KA-1]  Exp.  Coefficient  of  air. 

B4(p)=l/Tf4(p); 


%  Emissivity  of  titanium. 

E4(p)=abs((-5.0412e-008*T4(p)A2)+(2.4949e-004*T4(p))+0.051504); 

E5(p)=abs((-5.0412e-008*T5(p)A2)+(2.4949e-004*T5(p))+0.051504); 
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%  Solves  for  Rayleigh  and  Grashof  #'s. 
Gr4(p)=((g*B4(p)*(T4(p)-T5(p))*(LcA3))/(nu4(p)A2)); 

Ra4(p)=Gr4(p)*Pr4(p); 

ifRa4(p)=0 

Ra4(p)=0.001; 

else 

Ra4(p)=Ra4(p); 

end 

%  Nusselt  correlation  for  circular  annuli  heated  and  cooled  on  the 
%  vertical  curved  surface.  Developed  by  de  Vahl  Davis  and  Thomas. 
fpr4(p)=(  1  +((0 . 5/Pr4(p))A(9/ 1 6)))  A(- 1 6/9); 
Nu4(p)=0.364*((Ra4(p)*fpr4(p))A0.25)*((r3/r2)A0.5); 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  within  the  annuli. 
h4(p)=(Nu4(p)*ka4(p))/Lc; 

Rco2(p)=(A4+A5)/(h4(p)*A4*A5); 

%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  interior  of  outer  right  circular  cylinder 
%  of  finite  length  to  exterior  of  inner  right  circular  coaxial 
%  cylinder  [C-87], 

R=r3/r2; 

HH=Lc/r2; 

F45=(  1  /R)  *  ( 1  -((HHA2+RA2- 1  )/(4  *HH))- 1  /pi  *  (acos((HHA2-RA2+ 1  )/(HHA2+RA2- 1  ))-(((sq 
rt(((HHA2+RA2+l)A2)-(2*(RA2))))/(2*HH))*acos((HHA2-RA2+ 1  )/(R*(HHA2+RA2-l ))))-( 
((HHA2-RA2+ 1 )/ (2  *HH))  *  asin(  1  /R)))); 

%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

Rr2(p)=((l~E4(p))/(E4(p)*A4))+(l/(A4*F45))+((l-E5(p))/(E5(p)*A5)); 

Rrad2(p)=  1  /  ((sig/Rr2(p))  *  ((T  4(p)+T5(p))  *(T  4(p)A2+T  5  (p)A2»); 

%  Solves  for  total  resistance. 

Rt2(p)=((Rco2(p)*Rrad2(p))/(Rco2(p)+Rrad2(p))); 

%  Finite  difference  equation— explicit  method. 

T4(p+l)=2*(((H2(p)/Rc2(p))*(T3(p)-T4(p)))+((H2(p)/Rt2(p))*(T5(p)-T4(p)))+((T3(p)+ 

T4(p))/2))-T3(p+l); 
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if  T4(p+1)<300 

T4(p+1)=300; 

else 

T4(p+ 1  )=T4(p+ 1 ); 
end 


%%%  SOLVES  SECTION  4  [T5]  %%% 

%  [K]  Film  temperature  to  evaluate  properties  of  titanium. 

Tf5(pHT5(p)+T6(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  canister  wall. 

kc5(p)=abs((l  .70 1  le-008*Tf5(p)A3)-(4.33 16e-005*Tf5(p)A2)+(2.7279e-002*Tf5(p))+20. 
358); 

%  [J/kg*K]  Specific  Heat  of  canister  wall. 

cp5(p)=abs(((4.6480e-007)*Tf5(p)A3)-((1.4813e-003)*Tf5(p)A2)+(1.5387*Tf5(p))+139.9 

3); 

%  [mA2/s]  Thermal  diffusivity  of  canister  wall. 

alf5(p)=abs(((-4.532e-015)*Tf5(p)A3)+((1.7281e-01  l)*Tf5(p)A2)-(1.929e-008*Tf5(p))+l 
.3594e-005); 

%  [kg/mA3]  Density  of  canister  wall. 
rho5(p)=kc5(p)/(alf5(p)*cp5(p)); 

%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 
Rc3(p)=log(r2/rl)/(2*pi*kc5(p)*Lc); 

%  Finite  difference  equation-explicit  method. 

H3(p)=delt/(rho5(p)*cp5(p)*V3); 

TS(p+l)=2*(((H3(p)/Rt2(p))»(T4(p)-T5(p)))+((H3(p)/Rc3(p))*(T6(p)-T5(p)))+((T5(p)+ 

T6(p)V2))-T6(p); 


%%%  SOLVES  SECTION  5  [T6]  %%% 

%  [K]  Film  temperature  to  evaluate  air  properties  in  2nd  annulus. 
Tf6(p)=(T6(p)+T7(p)y2; 
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%  [W/m*K]  Thermal  conductivity. 

ka6(p)=abs((4.9198e-01  l*Tf6(p)A3)-(1.608e-007*Tf6(p)A2)+(2.014e-004*Tf6(p))-0.020 
796); 

%  [mA2/s]  Kinematic  viscosity. 

nu6(p)=abs(((7.786e-005*Tf6(p)A2)+(4.4053e-002*Tf6(p))-2.4972)*1.0e-006); 

%  Prandtl  #. 

Pr6(p)=abs((-2.5065e-011*Tf6(p)A3)+(8.4944e-008*Tf6(p)A2)-(1.0417e-004*Tf6(p))+0. 

74553); 

%  [KA-1]  Exp.  Coefficient  of  air. 

B6(p)=l/Tf6(p); 

%  Emissivity  of  titanium. 

E6(p)=abs((-5.0412e-008*T6(p)A2)+(2.4949e-004*T6(p))+0.051504); 

E7(p)=abs((-5.0412e-008*T7(p)A2)+(2.4949e-004*T7(p))+0.051504); 

%  Solves  for  Rayleigh  and  Grashof  #'s. 
Gr5(p)=(g*B6(p)*(T6(p)-T7(p))*(LcA3))/(nu6(p)A2); 

Ra5(p)=Gr5(p)*Pr6(p); 

ifRa5(p)=0 

Ra5(p)=0.001; 

else 

Ra5(p)=Ra5(p); 

end 

%  Nusselt  correlation  for  circular  annuli  heated  and  cooled  on  the 
%  vertical  curved  surface.  Developed  by  de  Vahl  Davis  and  Thomas. 
fpr5(p)=(l+((0.5/Pr6(p))A(9/ 1 6)))A(- 1 6/9); 
Nu5(p)=0.364*((Ra5(p)*fpr5(p))A0.25)*((rl/ro)A0.5); 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  within  the  annuli. 
h5(p)=Nu5(p)*ka6(p)/Lc; 

Rco3(p)=(A6+A7)/(h5(p)*A6*A7); 

%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  interior  of  outer  right  circular  cylinder 
%  of  finite  length  to  exterior  of  inner  right  circular  coaxial 
%  cylinder  [C-87], 

Rl=rl/ro; 
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HHl=Lc/ro; 

F67=(l/Rl)*(l-((HHlA2+RlA2-l)/(4*HHl))-l/pi*(acos((HHlA2-RlA2+l)/(HHlA2+RlA 
2-l))-(((sqrt(((HHlA2+RlA2+l)A2)-(2*(RlA2))))/(2*HHl))*acos((HHlA2-RlA2+l)/(Rl* 
(HH 1 A2+R 1 A2- 1  ))))-(((HH  1 A2-R 1 A2+ 1  )/(2  *HH  1 ))  *asin(  1  /R 1 )))); 

%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

Rr  3  (p)=((  1 -E6  (p))/(E6(p)  *  A6))+(  1  /( A7  *  F67))+((  1 -E7  (p))/(E7(p)  *  A7)) ; 

Rrad3(p)=  l/((sig^r3(p))*((T6(p)+T7(p))*(T6(p)A2+T7(p)A2))); 

%  Solves  for  total  resistance. 

Rt3  (p)=(Rco3  (p)*Rrad3  (p))/(Rco3  (p)+Rrad3(p)); 

%  Finite  difference  equation— explicit  method. 

T6(p+ 1  )=2  *(((H3  (p)/Rc3  (p))*(T  5  (p)-T  6(p)))+((H3  (p)/Rt3  (p))  *  (T7(p)-T  6(p)))+((T  5  (p)+ 
T  6(p))/2))-T5(p+ 1 ); 

if  T6(p+1)<300 

T6(p+1)=300; 

else 

T6(p+ 1  )=T  6(p+ 1 ); 
end 

%  [K]  Film  temperature  to  evaluate  properties  of  aluminum. 

Tf7(p)=(T7(p)+T8(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  missile. 
kc7(p)=237; 

%  [J/kg*K]  Specific  Heat  of  missile. 
cp7(p)=903; 

%  [mA2/s]  Thermal  diffusivity  of  missile. 
alf7(p)=97. 1  e-006; 

%  [kg/mA3]  Density  of  missile. 
rho7(p)=2702; 

%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 

Rc4(p)=log(ro/r)/(2  *pi  *  kc7(p)  *Lc); 
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%  Finite  difference  equation— explicit  method. 

H4(p)=(delt/(rho7(p)*cp7(p)*V4)); 

T7  (p+ 1  )=2  *  (((H4(p)/Rt3  (p))*  (T  6(p)-T7  (p)))+((H4(p)/Rc4(p))*  (T8(p)-T7  (p)))+((T7  (p)+ 
T8(p))/2))-T8(p); 

T8(p+l)=2*(((H4(p)/Rc4(p))*(T7(p)-T8(p)))+((T7(p)+T8(p))/2))-T7(p+l); 

if  T8(p+1)<300 

T8(p+1)=300; 

else 

T8(p+l)=T8(p+l); 

end 

end 

t=t/60; 

plot(t,T8);  grid; 
xlabel('Time  (min)'); 
ylabel('Temperature  (K)'); 
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APPENDIX  B.  CCL-FIRE  PROGRAM  CODE  (2B-CENTER) 

This  MATLAB  code  solve  for  radiation  resistance  using  the  Taylor  series  expansion  as 


described  in  the  radiation  section  of  Chapter  II.  This  code  is  set  up  for  the  center  CCL  and 

scenario  2. 

format  bank 
clear 


%  DIMENSIONS  FOR  CANISTERS  &  SURROUNDING  VLS  WALL. 


w=l; 

% 

Lw=l; 

% 

Lc=l; 

% 

Lm-1; 

% 

delx=0.0048; 

% 

r=0.00159; 

% 

ro=0.2586; 

% 

rl=0.2635; 

% 

r2=0.2667; 

% 

r3=0.3492; 

% 

r4=0.3571; 

% 

s=0.8382; 

% 

g=9.81; 

% 

sig=5.67e-008; 

% 

Al=w*Lw; 

% 

A2=A1; 

A3=2*pi*r4*Lc; 

% 

A4=2*pi*r3*Lc; 

% 

A5=2*pi*r2*Lc; 

% 

A6=2*pi*rl*Lc; 

% 

A7=2*pi*ro*Lc; 

% 

Vl=delx*w*Lw; 

V2=pi*Lc*(r4A2-r3 

A2); 

V3=pi*Lc*(r2A2-rl 

A2); 

V4=pi*roA2*Lm; 

[m]  Width  of  for  and  aft  bulkhead. 

[m]  Height  of  wall  surrounding  CCL  module, 
[m]  Height  of  inner/outer  canister  wall. 

[m]  Length  of  missile. 

[m]  Thickness  of  wall  surrounding  CCL  module, 
[m]  Near  center  of  missile  (1/16  in). 

[m]  Radius  of  missile. 

[m]  Inside  radius  of  inner  cylinder  wall. 

[m]  Outside  radius  of  inner  cylinder  wall. 

[m]  Inside  radius  of  outer  cylinder  wall. 

[m]  Outside  radius  of  outer  cylinder  wall. 

[m]  Distance  between  cylinder  CL. 

[m/sA2]  Gravitational  constant. 

[W/mA2*KA4]  Stefan-Boltzmann  constant. 
[mA2]Area  of  VLS  wall  perp  to  dir  of  heat  flux. 

[mA2]  Outer  surface  area  of  1  st  cylinder. 

[mA2]  Inner  surface  area  of  1st  cylinder. 

[mA2]  Outer  surface  area  of  2nd  cylinder. 

[mA2]  Inner  surface  area  of  2nd  cylinder. 

[mA2]  Surface  area  of  missile. 


%  [K]  INITIAL  AMBIENT  CONDITIONS  TO  START  ITERATION. 
T2(l)=300; 
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T3(l)=300; 

T4(l)=300; 

T5(l)=300; 

T6(l)=300; 

T7(l)=300; 

T8(l)=300; 


%  DEFINES  TIME  STEP  AND  NUMBER  OF  ITERATIONS 

S=60000;  %  Number  of  iterations. 

delt=0.36;  %  [s]  Time  step. 

%  BEGINS  LOOP  TO  SOLVE  FOR  ALL  TEMPS  IN  RESISTANCE  NETWORK. 

for  p=l:S 
pl=p 


%  [s]  Time. 
t(pMdelt*(p-l)); 

dT(p)=t(p)*  1.3333; 
dT  A(p)=t(p)  *0.0128; 

T 1  (p)=3 00+dT  (p); 

T 1  A(p)=  1 200+dT  A(p); 

if  Tl(p)<=1200 
Tl(p)=T!(p); 

Tl(p)-TlA(p); 

end 

t(p+l)=(delt*(p)); 
dT(p+l)=t(p+l)*  1.3333; 
dTA(p+l)=t(p+l)*0.0128; 

T 1  (p+ 1 )=3 00+dT  (p+ 1 ); 
TlA(p+l)=1200+dTA(p+l); 

if  Tl(p+I)<=1200 
Tl(p+l)=Tl(p+l); 
else 

Tl(p+l)=TlA(p+l); 

end 
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%%%  SOLVES  SECTION  1  [T2]  %%% 


%  [K]  Film  temperature  to  evaluate  properties  of  mild  steel 
Tfl(p)=(Tl(p)+T2(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  bulkhead. 

kwl(p)=abs((1.8869e-005*Tfl(p)A2)-(6.9896e-002*Tfl(p))+83.438); 

%  [J/kg*K]  Specific  Heat  of  bulkhead. 

cpl(p)~abs(((-2.2828e-006)*Tfl(p)A3)+((4.6812e-003)*Tfl(p)A2)-(2.3469*Tfl(p))+800.99); 

%  [mA2/s]  Thermal  diffusivity  of  bulkhead. 

alfl(p)=abs(((4.0118e-014)*Tfl(p)A3)-((7.9337e-01  l)*Tfl(p)A2)+(3.3733e-008*Tfl(p))+1.0388 
e-005); 

%  [kg/mA3]  Density  of  bulkhead, 
rho  1  (p)=kwl(p)/(cp  1  (p)*alfl  (p)); 

%  Solves  for  the  conductive  resistance  in  bulkhead. 

Rc  1  (p)=delx/  (kw  1  (p)  *  A2); 

%  [K]  Film  temperature  to  evaluate  properties  in  1st  air  space 
Tf2(p)=(T2(p)+T3(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  air. 

ka2(p)=abs((4.9198e-011*Tf2(p)A3)-(1.608e-007*Tf2(p)A2)+(2.014e-004*Tf2(p))-0.020796); 

%  [mA2/s]  Kinematic  viscosity  of  air. 

nu2(p)=abs(((7.786e-005*Tf2(p)A2)+(4.4053e-002*T£2(p))-2.4972)*1.0e-006); 

%  Prandtl  #. 

Pr2(p)=abs((-2.5065e-01 1  *Tf2(p)A3)+(8.4944e-008*Tf2(p)A2)-(l.  04 17e-004*Tf2(p))+0. 74553); 

%  [KA-1]  Exp.  Coefficient  of  air  . 

B2(p)=l/Tf2(p); 

%  Emissivity  of  mild  steel. 

E2(p)=abs((-9.8416e-008*T2(p)A2)+(2.6912e-004*T2(p))+0. 73699); 

%  Emissivity  of  titanium. 

E3  (p)=abs((-5 ,0412e-008*T3  (p)A2)+(2 . 4949e-004 *  T3  (p))+0 .051504); 
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%  Solves  for  Rayleigh  and  Grashof  #'s  WRT  bulkhead  inner  surface. 

Gr2(p)=((g*B2(p)  *  (T2(p)-T  £2(p))  *  (LwA3  ))/ (nu2(p)A2)); 

Ra2(p)=Gr2(p)*Pr2(p); 

%  Correlation  to  solve  for  the  average  Nusselt  #  over  a 
%  vertical  flat  plat.  Recommended  by  Churchill  and  Chu,  it  is 
%  good  over  entire  range  of  Rayleigh  #'s. 

Nu2(p)=(0.825+((0.387*(Ra2(p)A(l/6)))/((l+((0.492/Pr2(p))A(9/16)))A(8/27))))A2; 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  over  the  bulkhead  inner  surface. 
h2(p)=Nu2(p)*ka2(p)/Lw; 

Rco  1  a(p)=  1  /  (h2(p)  *  A2); 

%  Solves  for  Rayleigh  and  Grashof  #'s  WRT  cylinder  outer  surface. 
Gr3(p)=(g*B2(p)*(Tf2(p)-T3(p))*(LcA3))/(nu2(p)A2); 

Ra3(p)=Gr3(p)*Pr2(p); 

%  IF-THEN  LOOP  TO  DETERMINE  IF  THE  OUTER  WALL  OF  THE  CYLINDER 
%  CAN  BE  APPROXIMATED  AS  A  PLANE  WALL.  IF  A  PLANE  WALL  CHURCHILL 
%  AND  CHU’S  CORRELATION  FOR  A  VERTICAL  SURFACE  IS  USED.  IF 
%  THERE  ARE  APPRECIABLE  CURVATURE  EFFECTS  A  CORRELATION  DEVELOPED 
%  BY  LE  FEVRE  AND  EDE  BASED  ON  CYLINDER  HEIGHT. 

if  Gr3(p)==0 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)A2)/(5*(20+21*Pr2(p))))A0.25)+((4*Lc*(272+315*Pr2(p)))/(3 

5  *  2  *r4  *  (64+63  *Pr2(p)))); 

else 

if  ((2*r4)/Lc)>=35/(Gr3(p)Al/4) 

Nu3(p)=(0.825+((0.387*(Ra3(p)A(l/6)))/((l+((0.492/Pr2(p))A(9/16)))A(8/27))))A2; 

else 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)A2)/(5*(20+21*Pr2(p))))A0.25)+((4*Lc*(272+315*Pr(p)»/(35 

*2*r4*(64+63  *Pr2(p)))); 

end 

end 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  over  outer  surface  of  the  cylinder. 
h3  (p)=(Nu3  (p)  *ka2(p))/Lc; 

Rcolb(p)=l/(h3(p)*A3); 

%  Adds  convective  resistances  in  series. 

Rco  1  (p)=Rco  1  a(p)+Rco  1  b(p); 
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%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  infinite  plane  to  a  row  of  pipes  [C-6], 
D=(2*r4)/s; 

F23=l  -(( 1  -DA2)A0.5)+(D*atan(((  1  -DA2)/DA2)A0. 5)); 

%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

Rrl(p)=((l-E2(p))/(E2(p)*A2)Hl/(A2*F23))+((l-E3(p)V(E3(p)*A3)); 
q°  1  (p)=(sig*(T2(p)A4-T3(p)A4))/Rr  l(p); 


if  qol(p)==0 

Rradl(p)=Rrl(p)/(4*sig*T3(p)A3); 

else 

Rrad  1  (p)=  1  /((qo  1  (p)/ (T2(p)-T3  (p)))+((4  *  sig*  T3  (p)A3  )/Rr  1  (p))); 
end 

%  Solves  for  total  resistance,  [T3],  and  begins  next  iteration. 

Rt(p)=(Rco  1  (p)*Rrad  1  (p))/(Rrad  1  (p)+Rco  1  (p)); 

%  Finite  difference  equation— explicit  method. 

H 1  (p)=delt/(rho  l(p)*cpl(p)*Vl); 

T2(p+ 1  )=2  *(((H  1  (p)/Rc  1  (p))* (T 1  (p)-T2(p)))+((H  1  (p)/Rt(p))  *  (T3  (p)-T2(p)))+((T2(p)+T  1  (p))/2) 
)-Tl(p+l); 

if  T2(p+1)<300 

T2(p+1)=300; 

else 

T2(p+l)=T2(p+l); 

end 


%%%  SOLVES  SECTION  2  [T3]  %%% 

%  [K]  Film  temperature  to  evaluate  properties  of  titanium. 

Tf3(p)=(T3(p)+T4(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  canister  wall. 

kc3(p)=abs((1.7011e-008*TD(p)A3)-(4.3316e-005*TD(p)A2)+(2.7279e-002*Tf3(p))+20.358); 
%  [J/kg*K]  Specific  Heat  of  canister  wall. 

cp3(p)=abs(((4.6480e-007)*Tf3(p)A3)-((1.4813e-003)*Tf3(p)A2)+(1.5387*Tf3(p))+139.93); 
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%  [mA2/s]  Thermal  diffusivity  of  canister  wall. 

alf3(p)=abs(((-4.532e-015)*Tf3(p)A3)+((1.7281e-01  l)*Tf3(p)A2)-(1.929e-008*Tf3(p))+1.3594e- 
005); 

%  [kg/mA3]  Density  of  canister  wall. 
rho3  (p)=kc3  (p)/ (alf3  (p)  *  cp3  (p)); 

%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 

Rc2(p)=log(r4/ r3)/ (2  *  pi  *kc3  (p)  *Lc); 

%  Finite  difference  equation— explicit  method. 

H2(p)=delt/(rho3(p)*cp3(p)*V2); 

T3  (p+ 1  )=2  *(((H2(p)/Rt(p))  *  (T2(p)-T3  (p)))+((H2(p)/Rc2(p))  *  (T4(p)-T3  (p)))+((T3  (p)+T  4(p))/2) 
)-T4(p); 


%%%  SOLVES  SECTION  3  [T4]  %%% 

%  [K]  Film  temperature  to  evaluate  air  properties  in  1st  annulus. 

Tf4(p)=(T4(p)+T5(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  air. 

ka4(p):=abs((4.9198e-011*Tf4(p)A3)-(1.608e-007*Tf4(p)A2)+(2.014e-004*Tf4(p))-0. 020796); 

%  [mA2/s]  Kinematic  viscosity  of  air. 

nu4(p)=abs(((7.786e-005*Tf4(p)A2)+(4.4053e-002*Tf4(p))-2.4972)*1.0e-006); 

%  Prandtl  #. 

Pr4(p)=abs((-2.5065e-011*Tf4(p)A3)+(8.4944e-008*Tf4(p)A2)-(1.0417e-004*Tf4(p))+0.74553); 

%  [KA-1]  Exp.  Coefficient  of  air. 

B4(p)=l/Tf4(p); 

%  Emissivity  of  titanium. 

E4(p)=abs((-5.0412e-008*T4(p)A2)+(2.4949e-004*T4(p))+0.051504); 

E5(p)=abs((-5.0412e-008*T5(p)A2)+(2.4949e-004*T5(p))+0.051504); 

%  Solves  for  Rayleigh  and  Grashof  #'s. 

Gr4(p)=((g*B4(p)*(T4(p)-T5(p))*(LcA3))/(nu4(p)A2)); 

Ra4(p)=Gr4(p)*Pr4(p); 
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ifRa4(p)==0 

Ra4(p)=0.001; 

else 

Ra4(p)=Ra4(p); 

end 

%  Nusselt  correlation  for  curcular  annuli  heated  and  cooled  on  the 
%  vertical  curved  surface.  Developed  by  de  Vahl  Davis  and  Thomas. 
fpr4(p)=(  1  +((0. 5/Pr4(p))A(9/ 1 6)))A(- 16/9); 
Nu4(p)=0.364*((Ra4(p)*fpr4(p))A0.25)*((r3/r2)A0.5); 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  within  the  annuli. 
h4(p)=(Nu4(p)  *  ka4(p))/Lc; 

Rco2(p)=(  A4+ A5  )/(h4(p)  *  A4  *  A5); 

%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  interior  of  outer  right  circular  cylinder 
%  of  finite  length  to  exterior  of  inner  right  circular  coaxial 
%  cylinder  [C-87], 

R=r3/r2; 

HH=Lc/r2; 

F45  =(  1  /R)  *  ( 1  -((HHA2+RA2- 1 )/ (4  *HH))- 1  /pi  *  (acos((HHA2-RA2+ 1  )/(HFT 
HA2+RA2+l)A2)-(2*(RA2))))/(2*HH))*acos((HHA2-RA2+l)/(R*(HHA2+: 
+l)/(2*HH))*asin(l/R)»); 

%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

Rr2(p)=((l-E4(p))/(E4(p)*A4))+(l/(A4*F45))+((l-E5(p))/(E5(p)*A5)); 
q°2(p)=(sig*(T  4(p)A4-T5(p)A4))/Rr2(p); 

if  qo2(p)==0 

Rrad2(p)=Rr2(p)/(4*sig*T5(p)A3); 

else 

Rrad2(p)=  1  / ((qo2(p)/ (T  4(p)-T  5  (p)))+((4*  sig*T5  (p)A3  )/Rr2(p))); 
end 

%  Solves  for  total  resistance. 
Rt2(p)=((Rco2(p)*Rrad2(p))/(Rco2(p)+Rrad2(p))); 


'2+RA2- 1  ))-(((sqrt(((H 
lA2- 1  ))))-(((HHA2-RA2 


75 


%  Finite  difference  equation— explicit  method. 

T4(p+l)=2*(((H2(p)/RC2(p))*(T3(p)-T4(p)))+((H2(p)/Rt2(p))*(T5(p)-T4(p)))+((T3(p)+T4(p))/ 

2))-T3(p+l); 


if  T4(p+1)<300 

T4(p+1)=300; 

else 

T4(p+ 1  )=T  4(p+ 1 ); 
end 


%%%  SOLVES  SECTION  4  [T5]  %%% 

%  [K]  Film  temperature  to  evaluate  properties  of  titanium. 

Tf5(p)=(T5(p)+T6(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  canister  wall. 

kc5(p)=abs((1.7011e-008*Tf5(p)A3)-(4.3316e-005*Tf5(p)A2)+(2.7279e-002*Tf5(p))+20.358); 

%  [J/kg*K]  Specific  Heat  of  canister  wall. 

cp5(p)=abs(((4.6480e-007)*Tf5(p)A3)-((1.4813e-003)*Tf5(p)A2)+(1.5387*Tf5(p))+139.93); 

%  [mA2/s]  Thermal  diffusivity  of  canister  wall. 

alf5(p)=abs(((-4. 532e-0 1 5)*Tf5(p)A3)+((  1 . 728  le-0 1 1  )*Tf5(p)A2)-(  1 . 929e-008*Tf5(p))+l  .3  594e- 
005); 

%  [kg/mA3]  Density  of  canister  wall. 
rho5(p)=kc5(p)/(alf5(p)*cp5(p)); 

%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 
Rc3(p)=log(r2/rl)/(2*pi*kc5(p)*Lc); 

%  Finite  difference  equation— explicit  method. 

H3(p)=delt/(rho5(p)*cp5(p)*V3); 

T5  (p+ 1  )=2  *  (((H3  (p)/Rt2(p))  *  (T4(p)-T5(p)))+((H3  (p)/Rc3  (p))  *(T6(p)-T  5  (p)))+((T  5  (p)+T  6(p))/ 
2))-T6(p); 


%%%  SOLVES  SECTION  5  [T6]  %%% 

%  [K]  Film  temperature  to  evaluate  air  properties  in  2nd  annulus. 
TfB(p)=(T6(p)+T7(p))/2; 
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%  [W/m*K]  Thermal  conductivity. 

ka6(p)=abs((4.9198e-01  l*Tf6(p)A3)-(1.608e-007*Tf6(p)A2)+(2.014e-004*Tf6(p))-0. 020796); 

%  [mA2/s]  Kinematic  viscosity. 

nu6(p)=abs(((7.786e-005*Tf6(p)A2)+(4.4053e-002*Tf6(p))-2.4972)*1.0e-006); 

%  Prandtl  #. 

Pr6(p)=abs((-2.5065e-011*Tf6(p)A3)+(8.4944e-008*Tf6(p)A2)-(1.0417e-004*Tf6(p))+0. 74553); 

%  [KA-1]  Exp.  Coefficient  of  air. 

B6(p)=l/Tf6(p); 

%  Emissivity  of  titanium. 

E6(p)=abs((-5.0412e-008*T6(p)A2)+(2.4949e-004*T6(p))+0.051504); 

E7(p)=abs((-5 .041 2e-008*T7(p)A2)+(2.4949e-004*T7(p))+0.05 1 504); 

%  Solves  for  Rayleigh  and  Grashof  #'s. 

Gr5  (p)=(g*B6(p)  *  (T6(p)-T7  (p))*  (LcA3  ))/ (nu6(p)A2); 

Ra5(p)=Gr5(p)*Pr6(p); 


if  Ra5(p)==0 

Ra5(p)=0.001; 

else 

Ra5(p)=Ra5(p); 

end 

%  Nusselt  correlation  for  curcular  annuli  heated  and  cooled  on  the 
%  vertical  curved  surface.  Developed  by  de  Vahl  Davis  and  Thomas. 
fpr5(p)=(  1  +((0. 5/Pr6(p))A(9/l  6)))A(-1 6/9); 
Nu5(p)=0.364*((Ra5(p)*fpr5(p))A0.25)*((rl/ro)A0.5); 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  within  the  annuli. 
h5  (p)=Nu5(p)  *  ka6(p)/Lc; 

Rco3  (p)=(  A6+ A7)/ (h  5  (p)  *  A6  *  A7); 

%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  interior  of  outer  right  circular  cylinder 
%  of  finite  length  to  exterior  of  inner  right  circular  coaxial 
%  cylinder  [C-87], 

Rl=rl/ro; 

HHl=Lc/ro; 
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F67=(1/R1)*(1  -((HH 1 A2+R 1 A2- 1  )/(4  *HH  1 ))- 1  /pi  *  (acos((HH  1 A2-R 1 A2+ 1  )/(HH  1 A2+R1 A2- 1 ))-(( 
(sqrt(((HH  1 A2+R 1 A2+ 1  )A2)-(2  *(R  1 A2))))/ (2  *HH  1 ))  *  acos((HH  1 A2-R 1 A2+ 1  )/(R  1  *  (HH 1 A2+R 1 A2 
-  1))))-(((HH  1 A2-R1A2+ 1  )/(2*HHl  ))*asin(l/Rl )))); 

%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

Rr3(p)=((l-E6(p))/(E6(p)*A6)>t(l/(A7*F67))+((l.E7(p))/(E7(p)*A7)); 

qo3(p)=(sig*(T6(p)A4-T7(p)A4))/Rr3(p); 

if  qo3(p)==0 

Rrad3(p)=Rr3(p)/(4*sig*T7(p)A3); 

else 

Rrad3  (p)=  1  /((qo3  (p)/(T6(p)-T7  (p)))+((4  *  sig*  T7  (p)A3  )/Rr  3  (p))); 
end 

%  Solves  for  total  resistance. 

Rt3  (p)=(Rco3  (p)*Rrad3  (p))/(Rco3  (p)+Rrad3  (p)); 

%  Finite  difference  equation— explicit  method. 

T6(p+ 1  )=2  *  (((H3  (p)/Rc3  (p))  *(T5(p)-T  6(p)))+((H3  (p)/Rt3  (p))*  (T7  (p)-T  6(p)))+((T  5  (p)+T  6(p))/ 
2))-T5(p+l); 

ifT6(p+l)<300 

T6(p+1)=300; 

else 

T  6(p+ 1  )=T  6(p+ 1 ); 
end 

%  [K]  Film  temperature  to  evaluate  properties  of  aluminum. 

Tf7(p)=(T7(p)+T8(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  missile. 
kc7(p)=237; 

%  [J/kg*K]  Specific  Heat  of  missile. 
cp7(p)=903; 

%  [mA2/s]  Thermal  diffUsivity  of  missile. 
alf7(p)=97. 1  e-006; 

%  [kg/mA3]  Density  of  missile. 
rho7(p)=2702; 
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%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 
Rc4(p)=log(ro/r)/(2*pi*kc7(p)*Lc); 

%  Finite  difference  equation— explicit  method. 

H4(p)=(delt/(rho7(p)*cp7(p)*V4)); 

T7(p+l)=2*(((H4(p)/Rt3(p»*(T6(p)-T7(p)))+((H4(p)/Rc4(p))*(T8(p)-T7(p)))+((T7(p)+T8(p))/ 

2))-T8(p); 

T8(p+l)=2*(((H4(p)/Rc4(p))*(T7(p)-T8(p)))+((T7(p)+T8(p))/2))-T7(p+l); 

if  T8(p+1)<300 

T8(p+1)=300; 

else 

T8(p+ 1  )=T8(p+ 1 ); 

end 

end 

t=t/60; 

plot(t,T8);  grid; 
xlabel('Time  (min)'); 
ylabel('Temperature  (K)'); 
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APPENDIX  C.  CCL-FIRE  PROGRAM  CODE  (2A-CORNER) 

This  MATLAB  code  solves  for  the  radiation  resistance  explicitly  as  described  in 

the  radiation  section  of  Chapter  II.  This  code  is  set  up  for  scenario  2  and  the  corner  CCL 

in  which  the  adjoining  bulkhead  is  assumed  to  be  a  reradiating  surface. 

format  bank 
clear 


%  DIMENSIONS  FOR  CANISTERS  &  SURROUNDING  VLS  WALL. 


w=l; 

% 

wl=l; 

% 

Lw=l; 

% 

Lc=l; 

% 

Lm=l; 

% 

delx=0.0048; 

% 

r=0. 00159; 

% 

ro=0.2586; 

% 

rl=0.2635; 

% 

r2=0.2667; 

% 

r3=0.3492; 

% 

r4=0.3571; 

% 

s=0.8382; 

% 

g=9.81; 

% 

sig=5.67e-008; 

% 

AR=wl*Lw; 

% 

A2=w*Lw; 

% 

A3=2*pi*r4*Lc; 

% 

A4=2*pi*r3*Lc; 

% 

A5=2*pi*r2*Lc; 

% 

A6=:2*pi*rl*Lc; 

% 

A7=2*pi*ro*Lc; 

% 

Vl=delx*w*Lw; 

V2=pi*Lc*(r4A2-r3 

A2); 

V3=pi*Lc*(r2A2-rl 

A2); 

V4=pi*roA2*Lm; 

[m]  Width  of  for  and  aft  bulkhead. 

[m]  Width  of  port  and  stbd  bulkhead. 

[m]  Height  of  wall  surrounding  CCL  module. 

[m]  Height  of  inner/outer  canister  wall. 

[m]  Length  of  missile. 

[m]  Thickness  of  wall  surrounding  CCL  module, 
[m]  Near  center  of  missile  (1/16  in). 

[m]  Radius  of  missile. 

[m]  Inside  radius  of  inner  cylinder  wall. 

[m]  Outside  radius  of  inner  cylinder  wall. 

[m]  Inside  radius  of  outer  cylinder  wall. 

[m]  Outside  radius  of  outer  cylinder  wall. 

[m]  Distance  between  cylinder  CL. 

[m/sA2]  Gravitational  constant. 

[W/mA2*KA4]  Stefan-Boltzman  constant. 

[mA2]  Area  of  reradiative  bulkhead 

[mA2]  Area  of  VLS  wall  perp  to  dir  of  heat  flux. 

[mA2]  Outer  surface  area  of  1st  cylinder. 

[mA2]  Inner  surface  area  of  1st  cylinder. 

[mA2]  Outer  surface  area  of  2nd  cylinder. 

[mA2]  Inner  surface  area  of  2nd  cylinder. 

[mA2]  Surface  area  of  missile. 
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%  [K]  INITIAL  AMBIENT  CONDITIONS  TO  START  ITERATION. 


T2(l)=300; 

T3(l)=300; 

T4(l)=300; 

T5(l)=300; 

T6(l)-300; 

T7(l)=300; 

T8(l)=300; 

%  DEFINES  TIME  STEP  AND  NUMBER  OF  ITERATIONS 

S=45000;  %  Number  of  iterations. 

delt=0.36;  %  [s]  Time  step. 


%  BEGINS  LOOP  TO  SOLVE  FOR  ALL  TEMPS  IN  RESISTANCE  NETWORK. 


for  p=l:S 
pl=p 


%  [s]  Time. 

t(p)=(delt*(p-l)); 

%  Temperature  rise  from  ambient  to  1325  K  over  a  5  minute  period. 
%  This  temp  is  then  maintained  for  the  duration  of  the  fire. 
dT(p)=t(p)*  1.3333; 
dTA(p)=t(p)*0.0128; 

T 1  (p)=3 00+dT  (p) ; 

T 1  A(p)=  1 200+dT  A(p); 

if  Tl(p)<=1200 
Tl(p)=Tl(p); 

Tl(p)=TlA(p), 

end 

t(p+l)=(delt*(p)); 
dT(p+l)=t(p+l)*  1.3333; 
dTA(p+l)=t(p+l)*0.0128; 

T 1  (p+ 1  )=3  00+dT  (p+ 1 ); 

T 1  A(p+ 1  )=1 200+dTA(p+ 1 ); 
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if  Tl(p+l)<=1200 
Tl(p+l)=Tl(p+l); 
else 

T 1  (p+ 1  )=T  1  A(p+ 1 ); 
end 


%%%  SOLVES  SECTION  1  [T2]  %%% 

%  [K]  Film  temperature  to  evaluate  properties  of  mild  steel 
Ttl(pHTl(p)+T2(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  bulkhead. 
kwl(p)=abs((1.8869e-005*Tfl(p)A2)-(6.9896e-002*Tfl(p))+83.438); 

%  [J/kg*K]  Specific  Heat  of  bulkhead. 

cpl(p)=abs(((-2.2828e-006)*Tfl(p)A3)+((4.6812e-003)*Tfl(p)A2)-(2.3469*Tfl(p))+800. 

99); 

%  [mA2/s]  Thermal  diffusivity  of  bulkhead. 

alfl(p)=abs(((4.0118e-014)*Tfl(p)A3)-((7.9337e-011)*Tfl(p)A2)+(3.3733e-008*Tfl(p)) 

+1.0388e-005); 

%  [kg/mA3]  Density  of  bulkhead, 
rho  1  (p)=kw  1  (p)/(cp  1  (p)  *  alfl  (p)); 

%  Solves  for  the  conductive  resistance  in  bulkhead. 

Rc  1  (p)=delx/  (kw  1  (p)  *  A2) ; 

%  [K]  Film  temperature  to  evaluate  properties  in  1st  air  space. 

Tf2(p)=(T2(p)+T3(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  air. 

ka2(p)=abs((4.9198e-01 1  *Tf2(p)A3)-(l  ,608e-007*Tf2(p)A2)+(2.014e-004*Tf2(p))-0.020 
796); 

%  [mA2/s]  Kinematic  viscosity  of  air. 

nu2(p)=abs(((7.786e-005*Tf2(p)A2)+(4.4053e-002*Tf2(p))-2.4972)*1.0e-006); 

%  Prandtl  #. 

Pr2(p)=abs((-2 .5065 e-0 1 1  *  T  f2(p)A3 )+(8 . 4944e-008 *Tf2(p)A2)-(  1 . 04 1 7  e-004 *Tf2(p))+0 . 
74553); 
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%  [K A- 1  ]  Exp.  Coefficient  of  air  . 

B2(p)=l/Tf2(p); 

%  Emissivity  of  mild  steel. 

E2(p)=abs((-9. 841 6e-008 *T2(p)A2)+(2 .691 2e-004*T2(p))+0. 73 699); 

%  Emissivity  of  titanium. 

E3(p)=abs((-5.0412e-008*T3(p)A2)+(2.4949e-004*T3(p))+0.051504); 

%  Solves  for  Rayleigh  and  Grashof  #'s  WRT  bulkhead  inner  surface. 
Gr2(p)=((g*B2(p)*(T2(p)-Tf2(p))*(LwA3))/(nu2(p)A2)); 

Ra2(p)=Gr2(p)*Pr2(p); 

%  Correlation  to  solve  for  the  average  Nusselt  #  over  a 
%  vertical  flat  plat.  Recommended  by  Churchill  and  Chu,  it  is 
%  good  over  entire  range  of  Rayleigh  #'s. 

Nu2(p)=(0  .  825+((0 . 3  8 7 * (Ra2(p) A(  1 /6)))/(( 1  +((0 . 492/Pr2(p))A(9/ 1 6)))A(8/27))))A2; 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  over  the  bulkhead  inner  surface. 
h2(p)=Nu2(p)*ka2(p)/Lw; 

Rco  1  a(p)=  1  /(h2(p)  *  A2); 

%  Solves  for  Rayleigh  and  Grashof  #'s  WRT  cylinder  outer  surface. 

Gr3  (p)=(g*B2(p)  *  (Tf2(p)-T3  (p))  *  (LcA3  ))/ (nu2(p)  A2); 

Ra3(p)=Gr3(p)*Pr2(p); 

%  IF-THEN  LOOP  TO  DETERMINE  IF  THE  OUTER  WALL  OF  THE  CYLINDER 
%  CAN  BE  APPROXIMATED  AS  A  PLANE  WALL.  IF  A  PLANE  WALL 
%  CHURCHILL 

%  AND  CHU’S  CORRELATION  FOR  A  VERTICAL  SURFACE  IS  USED.  IF 
%  THERE  ARE  APPRECIABLE  CURVATURE  EFFECTS  A  CORRELATION 
%  DEVELOPED 

%  BY  LE  FEVRE  AND  EDE  BASED  ON  CYLINDER  HEIGHT. 


if  Gr3(p)= 0 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)A2)/(5*(20+21*Pr2(p))))A0.25)+((4*Lc*(272+315*Pr2 

(p)))/ (3  5  *  2  *  r4 *  (64+63  *Pr2(p)))) ; 

else 

if  ((2*r4)/Lc)>=3  5/(Gr3  (p)Al/4) 

Nu3(p)=(0.825+((0.387*(Ra3(p)A(l/6)))/((l+((0.492/Pr2(p))A(9/16)))A(8/27))))A2; 

else 
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Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)A2)/(5*(20+21*Pr2(p))))A0.25)+((4*Lc*(272+3 15*Pr2 

(p)))/(35*2*r4*(64+63*Pr2(p))»; 

end 

end 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  over  outer  surface  of  the  cylinder. 
h3  (p)=(Nu3  (p)  *ka2(p))/Lc; 

Rcolb(p)=l/(h3(p)*A3); 

%  Adds  convective  resistances  in  series. 

Rco  1  (p)=Rco  1  a(p)+Rco  1  b(p); 

%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  infinite  plane  to  a  row  of  pipes  [C-6], 

D=(2*r4)/s; 

F23 = 1  -(( 1  -D  A2)A0 . 5  )+(D  *  atan(((  1  -D  A2)/D  A2)A0. 5 )) ; 

FR3=F23; 

F2R=0 . 5  *  ( 1  +(w  1  /w)-sqrt(  1  +((w  1  /w)  A2))); 

%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

J1=1/(A2*F2R); 

J2=1/(AR*FR3); 

JT=(J1+J2)A(-1); 

Rr  1  (p)=((  1  -E2(p))/ (E2(p)*  A2))+(  1  /  ((A2*F23  )+JT))+((  1 -E3  (p))/ (E3  (p)*  A3  )); 

Rrad  1  (p)=  1  /  ((sig/Rr  1  (p))  *  ((T2(p)+T3  (p))  *  (T2(p)  A2+T3  (p)A2))); 

%  Solves  for  total  resistance,  [T3],  and  begins  next  iteration. 

Rt(p)=(Rco  1  (p)*Rrad  1  (p))/(Rrad  1  (p)+Rco  1  (p»; 

%  Finite  difference  equation— explicit  method. 

H 1  (p)=delt/(rho  l(p)*cpl(p)*Vl); 

T2(p+ 1  )=2  *(((H1  (p)/Rc  1  (p))*(T  1  (p)-T2(p)))+((H  1  (p)/Rt(p))*(T3  (p)-T2(p)))+((T2(p)+T 
l(p))/2))-Tl(p+l); 

if  T2(p+1)<300 

T2(p+1)=300; 

else 

T2(p+ 1  )=T2(p+ 1 ); 
end 
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%%%  SOLVES  SECTION  2  [T3]  %%% 


%  [K]  Film  temperature  to  evaluate  properties  of  titanium. 

Tf3(p)=(T3(p)+T4(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  canister  wall. 

kc3(p)=abs((1.701  le-008*Tf3(p)A3)-(4.3316e-005*Tf3(p)A2)+(2.7279e-002*TD(p))+20 
358); 

%  [J/kg*K]  Specific  Heat  of  canister  wall. 

cp3(p)=abs(((4.6480e-007)*Tf3(p)A3)-((1.4813e-003)*Tf3(p)A2)+(1.5387*Tf3(p))+139.9 

3); 

%  [mA2/s]  Thermal  diffusivity  of  canister  wall. 

alf3(p)=:abs(((-4.532e-015)*Tf3(p)A3)+((1.7281e-011)*Tf3(p)A2)-(L929e-008*Tf3(p))+l 

,3594e-005); 

%  [kg/mA3]  Density  of  canister  wall. 
rho3  (p)=kc3  (p)/  (alf3  (p)  *  cp3  (p)) ; 

%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 
Rc2(p)=log(r4/r3)/(2*pi*kc3(p)*Lc); 

%  Finite  difference  equation— explicit  method. 

H2(p)=delt/ (rho3  (p)  *cp3  (p)  *  V2) ; 

T3(p+l)=2*(((H2(p)/Rt(p))*(T2(p)-T3(p)))+((H2(p)/Rc2(p))*(T4(p)-T3(p)))+((T3(p)+T 

4(p))/2))-T4(p); 


%%%  SOLVES  SECTION  3  [T4]  %%% 

%  [K]  Film  temperature  to  evaluate  air  properties  in  1st  annulus. 

T  f4(p)=(T  4(p)+T  5  (p»/2; 

%  [W/m*K]  Thermal  conductivity  of  air. 

ka4(p)=abs((4.9 1 98e-0 1 1  *Tf4(p)A3)-(1.608e-007*Tf4(p)A2)+(2.014e-004*Tf4(p))-0  020 
796); 

%  [mA2/s]  Kinematic  viscosity  of  air. 

nu4(p)=abs(((7.786e-005  *Tf4(p)A2)+(4.4053e-002*Tf4(p))-2.4972)*  1 .0e-006); 
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%  Prandtl  #. 

Pr4(p)=abs((-2 . 5065 e-0 1 1  *Tf4(p)  A3 )+(8 . 4944e-008 *  Tf4(p)A2)-(  1 ,0417e-004*Tf4(p))+0. 
74553); 

%  [KA-1]  Exp.  Coefficient  of  air. 

B4(p)=l/Tf4(p); 

%  Emissivity  of  titanium. 

E4(p)=abs((-5.04 1 2e-008*T4(p)A2)+(2.4949e-004*T4(p))+0.05 1 504); 
E5(p)=abs((-5.0412e-008*T5(p)A2)+(2.4949e-004*T5(p))+0.051504); 

%  Solves  for  Rayleigh  and  Grashof  #'s. 
Gr4(p)=((g*B4(p)*(T4(p)-T5(p))*(LcA3))/(nu4(p)A2)); 

Ra4(p)=Gr4(p)*Pr4(p); 

ifRa4(p)==0 

Ra4(p)=0.001; 

else 

Ra4(p)=Ra4(p); 

end 

%  Nusselt  correlation  for  circular  annuli  heated  and  cooled  on  the 
%  vertical  curved  surface.  Developed  by  de  Vahl  Davis  and  Thomas. 
fpr4(p)=(  1  +((0 . 5/Pr4(p))A(9/ 1 6)))  A(- 16/9); 
Nu4(p)=0.364*((Ra4(p)*fpr4(p))A0.25)*((r3/r2)A0.5); 

%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  within  the  annuli. 
h4(p)=(Nu4(p)*ka4(p))/Lc; 

Rco2(p)=(A4+A5)/(h4(p)*A4*A5); 

%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  interior  of  outer  right  circular  cylinder 
%  of  finite  length  to  exterior  of  inner  right  circular  coaxial 
%  cylinder  [C-87], 

R=r3/r2; 

HH=Lc/r2; 

F45=(l/R)*(l  -((HHA2+RA2- 1  )/(4*HH))-  l/pi*(acos((HHA2-RA2+  1)/(HHA2+RA2- 1  ))-(((sq 
rt(((HHA2+RA2+ 1  )A2)-(2*(RA2))))/(2*HH))*acos((HHA2-RA2+l)/(R*(HHA2+RA2-l))))-( 
((HHA2-RA2+ 1 )/ (2  *HH))  *  asin(  1  /R)))); 
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%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

Rr2(p)=((  1  -E4(p))/(E4(p)  *  A4))+(  1  /(A4  *F4  5))+((  1 -E5  (p))/ (E5  (p)  *  A5 )); 

Rrad2(p)=  l/((sig^r2(p))*((T4(p)+T5(p))*(T4(p)^2+T5(Pr2))); 

%  Solves  for  total  resistance. 

Rt2(p)=((Rco2(p)*Rrad2(p))/(Rco2(p)+Rrad2(p))); 

%  Finite  difference  equation--explicit  method. 

T4(p+ 1  )=2  *(((H2(p)/Rc2(p))*(T3  (p)-T4(p)))+((H2(p)/Rt2(p))  *(T5(p)-T4(p)))+((T3  (p)+ 
T  4(p))/2))-T3  (p+ 1 ); 

if  T4(p+1)<300 

T4(p+1)=300; 

else 

T  4(p+ 1  )=T  4(p+ 1 ); 
end 


%%%  SOLVES  SECTION  4  [T5]  %%% 

%  [K]  Film  temperature  to  evaluate  properties  of  titanium. 

Tf5(p)=(T5(p)+T6(p)y2; 

%  [W/m*K]  Thermal  conductivity  of  canister  wall. 

kc5(p)=abs((  1.701  le-008*Tf5(p)A3)-(4.3316e-005*Tf5(p)A2)+(2.7279e-002*Tf5(p))+20. 
358); 

%  [J/kg*K]  Specific  Heat  of  canister  wall. 

cp5(p)=abs(((4.6480e-007)*Tf5(p)A3)-((1.4813e-003)*Tf5(p)A2)+(1.5387*Tf5(p))+139.9 

3); 

%  [mA2/s]  Thermal  diffusivity  of  canister  wall. 

alf5(p)=abs(((-4.532e-015)*Tf5(p)A3)+((1.7281e-011)*Tf5(p)A2)-(1.929e-008*Tf5(p))+l 

,3594e-005); 

%  [kg/mA3]  Density  of  canister  wall. 
rho5(p)=kc5(p)/(alf5(p)*cp5(p)); 

%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 
Rc3(p)=log(r2/rl)/(2*pi*kc5(p)*Lc); 
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%  Finite  difference  equation— explicit  method. 

H3(p)=delt/(rho5(p)*cp5(p)*V3); 

T5(p+l)=2*(((H3(p)/Rt2(p))*(T4(p)-T5(p)))+((H3(p)/Rc3(p))*(T6(p)-T5(p)))+((T5(p)+ 

T6(p))/2»-T6(p); 


%%%  SOLVES  SECTION  5  [T6]  %%% 

%  [K]  Film  temperature  to  evaluate  air  properties  in  2nd  annulus. 
TfS(p)=(T6(p)+T7(p))/2; 

%  [W/m*K]  Thermal  conductivity. 

ka6(p)=abs((4.9198e-011*Tf6(p)A3)-(1.608e-007*Tf6(p)A2)+(2.014e-004*Tf6(p))-0.020 

796); 

%  [mA2/s]  Kinematic  viscosity. 

nu6(p)=abs(((7.786e-005*Tf6(p)A2)+(4.4053e-002*Tf6(p))-2.4972)*1.0e-006); 

%  Prandtl  #. 

Pr6(p)=abs((-2 . 5 06 5 e-0 1 1  *Tf6(p)A3)+(8.4944e-008*Tf6(p)A2)-(1.0417e-004*Tf6(p))+0. 
74553); 

%  [KA-1]  Exp.  Coefficient  of  air. 

B6(p)=l/Tf6(p); 

%  Emissivity  of  titanium. 

E6(p)=abs((-5 . 04 1 2e-008  *T6(p)A2)+(2 ,4949e-004 *  T6(p))+0 . 05 1 5 04); 

E7(p)=abs((-5 .04 1 2e-008*T7(p)A2)+(2.4949e-004*T7(p))+0.05 1 504); 

%  Solves  for  Rayleigh  and  Grashof  #'s. 
Gr5(p)=(g*B6(p)*(T6(p)-T7(p))*(LcA3))/(nu6(p)A2); 

Ra5(p)=Gr5(p)*Pr6(p); 

ifRa5(p)==0 

Ra5(p)=0.001; 

else 

Ra5(p)=Ra5(p); 

end 

%  Nusselt  correlation  for  circular  annuli  heated  and  cooled  on  the 
%  vertical  curved  surface.  Developed  by  de  Vahl  Davis  and  Thomas. 
fpr5(p)=(  1  +((0. 5/Pr6(p))A(9/ 1 6)))A(- 1 6/9); 
Nu5(p)=0.364*((Ra5(p)*fpr5(p))A0.25)*((rl/ro)A0.5); 
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%  Solves  for  convective  resistance  using  the  average  convection 
%  heat  transfer  coefficient  within  the  annuli. 
h5(p)=Nu5(p)*ka6(p)/Lc; 

Rco3(p)=(A6+A7)/(h5(p)*A6*A7); 

%  Calculates  view  factor.  Equation  from  the  catalog  of  radiation 
%  configuration  factors,  interior  of  outer  right  circular  cylinder 
%  of  finite  length  to  exterior  of  inner  right  circular  coaxial 
%  cylinder  [C-87], 

Rl=rl/ro; 

HHl=Lc/ro; 

F67=(l/Rl)*(l-((HHlA2+RlA2-l)/(4*HHl))-l/pi*(acos((HHlA2-RlA2+l)/(HHlA2+RlA 
2- 1  ))-(((sqrt(((HH  1 A2+R1 A2+ 1  )A2)-(2*  (R 1  A2))))/(2*HH  1  ))*acos((HH  1 A2-R1 A2+ 1  )/(R  1  * 
(HH 1 A2+R 1 A2- 1  ))))-(((HH  1 A2-R 1 A2+ 1 )/ (2  *HH  1 ))  *asin(  1  /R 1 )))) ; 

%  Solves  the  radiative  resistance  network  for  two  surf  enclosure, 

%  and  linearizes  the  resistance  to  be  used  ICW  the  convective 
%  resistance. 

Rr  3  (p)=((  1  -E6(p))/ (E6(p)  *  A6))+(  1  / ( A7  *F67))+((  1  -E7(p))/ (E7  (p)  *  A7)) ; 

Rrad3(p)=  l/((sig^r3fp))*((T6(p)+T7(p))*(T6(p)A2+T7(p)A2))); 

%  Solves  for  total  resistance. 

Rt3  (p)=(Rco3  (p)  *Rrad3  (p))/(Rco3  (p)+Rrad3  (p)) ; 

%  Finite  difference  equation— explicit  method. 

T  6(p+ 1  )=2  *(((H3  (p)/Rc3  (p))*  (T5  (p)-T6(p)))+((H3  (p)/Rt3  (p))  *  (T7  (p)-T  6(p)))+((T  5  (p)+ 
T  6(p))/2))-T5(p+ 1 ); 

if  T6(p+1)<300 

T6(p+1)=300; 

else 

T6(p+l)=T6(p+l); 

end 

%  [K]  Film  temperature  to  evaluate  properties  of  aluminum. 

Tf7(p)=(T7(p)+T8(p))/2; 

%  [W/m*K]  Thermal  conductivity  of  missile. 
kc7(p)=237; 

%  [J/kg*K]  Specific  Heat  of  missile. 
cp7(p)=903; 
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%  [mA2/s]  Thermal  diffusivity  of  missile. 
alf7(p)=97. 1  e-006; 


%  [kg/mA3]  Density  of  missile. 
rho7(p)=2702; 

%  Solves  for  the  conductive  resistance  through  the  cylinder  wall. 

Rc4(p)=log(ro/r)/(2  *  pi  *kc7(p)  *Lc); 

%  Finite  difference  equation-explicit  method. 

H4(p)=(delt/(rho7  (p)  *  cp7  (p)  *  V4)) ; 

T7(p+ 1  )=2  *  (((H4(p)/Rt3  (p))  *(T  6(p)-T7  (p)))+((H4(p)/Rc4(p))*  (T8  (p)-T7  (p)))+((T7  (p)+ 
T8(p))/2»-T8(p); 

T8(p+l)=2*(((H4(p)/Rc4(p))*(T7(p)-T8(p)))+((T7(p)+T8(p))/2))-T7(p+l); 

if  T8(p+1)<300 

T8(p+1)=300; 

else 

T8(p+l)=T8(p+l); 

end 

end 

t=t/60; 

plot(t,T8);  grid; 
xlabel('Time  (min)'); 
ylabel('Temperature  (K)'); 
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APPENDIX  D.  TIME-TEMP  PROFILE  AT  EACH  NODE 

The  following  figures  show  the  time-temperature  profile  for  each  node  in  the 
resistance  network.  These  four  figures  reflect  Scenario  1 .  Figures  24  and  25  reflect  the 
CCL-FIRE  program  code  which  solves  for  surface  radiative  resistance  explicitly.  Figures 
26  and  27  reflect  the  CCL-FIRE  program  code  which  solves  for  the  surface  radiative 
resistance  with  a  Taylor  series  expansion. 


Figure  24.  Scenario  1 A  -  Node  Time-Temperature  Profiles. 
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Figure  25.  Scenario  1A  -  Node  Time-Temperature  Profiles. 
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